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Abstract 

It was recently realized that quenched disorder may enhance the reliability of topological qubits 
by reducing the mobility of anyons at zero temperature. Here we compute storage times with and 
without disorder for quantum chains with unpaired Majorana fermions — the simplest toy model 
of a quantum memory. Disorder takes the form of a random site-dependent chemical potential. 
The corresponding one-particle problem is a one-dimensional Anderson model with disorder in the 
hopping amplitudes. We focus on the zero-temperature storage of a qubit encoded in the ground 
state of the Majorana chain. Storage and retrieval are modeled by a unitary evolution under the 
memory Hamiltonian with an unknown weak perturbation followed by an error-correction step. 
Assuming dynamical localization of the one-particle problem, we show that the storage time grows 
exponentially with the system size. We give supporting evidence for the required localization 
property by estimating Lyapunov exponents of the one-particle eigenfunctions. We also simulate 
the storage process for chains with a few hundred sites. Our numerical results indicate that in the 
absence of disorder, the storage time grows only as a logarithm of the system size. We provide 
numerical evidence for the beneficial effect of disorder on storage times and show that suitably 
chosen pseudorandom potentials can outperform random ones. 
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1 Introduction 



Topological features of physical systems are - by definition - robust against local imperfections. 
As a consequence, they are the natural object of study in the characterization and classification 
of different phases of matter. A major challenge is to identify and compute topological invariants for 
various systems and to connect these to experimentally accessible quantities. Examples include Chern 
numbers, the topological entanglement entropy [41, 43], or the classification of free fermions in terms 
of symmetry classes. 

In contrast, topological quantum computing attempts to make operational use of robust features of 
topologically ordered systems. A basic first goal is the reliable storage and readout of quantum infor- 
mation. Proposed toy models for topological quantum memories include Kitaev's Majorana chains [40] 
and the toric code [39]. These systems exhibit a 'topological' ground space degeneracy which persists 
in the presence of local perturbations. In both cases, the ground space constitutes a quantum error- 
correcting code. When superselection rules are taken into account (in the case of Majorana chains), 
these codes have a macroscopic distance which scales linearly with the system size. 

Unfortunately, these remarkable ground space properties are insufficient to protect the encoded 
information in the presence of generic local perturbations and/or thermal fiuctuations: localized exci- 
tations (or anyons) can be created at constant energy cost and once created, can propagate freely to 
cause logical errors. This fact was recognized early on [21] and has since been elaborated by several 
authors by giving explicit examples (see e.g., [48, 51]) and general impossibility proofs for storage 
at non-zero temperature. These range from the investigation of the robustness of topological order 
using finite-temperature generalizations of the topological entanglement entropy [13] to bounds on 
operational quantities such as the quantum memory relaxation rate [14] . For the 2D toric code, Alicki 
et al. showed that a weak coupling to a Markovian environment implies a constant lower bound on 
the relaxation rate of any logical state independent of the system size [3]. 

Several approaches for combating decoherence in quantum memories have been considered in the 
literature. Arguably the most intriguing of these is the idea of self-correcting quantum memories. This 
completely circumvents the decoherence mechanism mentioned above: it is based on a Hamiltonian 
which has no point-like excitations carrying a topological charge (such as anyons). More specifically, 
it assumes a Hamiltonian whose ground space is a quantum code with macroscopic distance and the 
property that high-weight errors correspond to larger energies. As a result, encoded states are sepa- 
rated by (macroscopic) energy barriers, a necessary condition for self-correction as formalized in [10] 
(see also [37]). A 4D generalization of the toric code is known [21] to provide a stable quantum memory 
in this fashion: it was shown in [4] that this memory is thermally stable in the sense that relaxation 
times increase exponentially with system size. Unfortunately, this model cannot be embedded in a 
locality-preserving way in M^. The existence of self-correcting memories in 3D remains under investi- 
gation and is of fundamental theoretical interest. On the other hand, simpler geometries, such as 2D 
(subspace) stabilizer codes cannot provide this kind of protection, as shown in [10]. 

A different line of research attempts to specifically address the creation and movement of anyons 
responsible for decoherence in topologically ordered systems. Importantly, while these excitations 
are thermally suppressed at low temperatures, their propagation impacts the coherence of quantum 
memories even in the case of zero temperature. Limiting this movement is therefore at the core of 
these methods. Each of these methods has its own merits and drawbacks. 

An example is active error correction which employs frequent rounds of syndrome measurements 
revealing the locations of anyons at any given time step. The syndrome information can then be used 
to undo the effect of errors caused by the propagation of anyons. It has been shown to permit reliable 
storage of quantum information for arbitrarily long times in e.g., the 2D toric code [21], with errors 
modeled as depolarizing noise. However, this approach requires significant experimental capabilities: 
the syndrome measurements need to be sufficiently accurate and sufficiently fast in order for the error 
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correction to succeed [21]. 

Given the stringent experimental requirements for active error-correction, it is natural to consider 
other mechanisms for reducing the propagation of anyons in topological quantum memories. One 
such approach is the introduction of long-range interactions between anyons. A concrete model where 
such an interaction is mediated by a bosonic background field was proposed by Hamma at al [28]. 
An alternative approach in which creation of anyons from the vacuum is penalized by long-range 
repulsive interactions was studied by Chesi et al [15]. A promising 3D stabilizer model was proposed 
more recently by Haah [27]. In this model anyons are bound together by energy barriers growing 
logarithmically with the lattice size [7] even though the underlying Hamiltonian includes only short- 
range interactions with constant strength. 

Quenched disorder in topological quantum memories 

Here we focus on another potentially more realistic approach towards limiting the detrimental effects 
of propagating anyons pioneered by Wootton and Pachos [56], and, independently, by Stark et al [51]. 
It relies on the observation that the mobility of anyons can be strongly reduced in the presence of 
quenched disorder. Such disorder may account for natural imperfections in physical realizations, or 
may be artificially engineered by tuning the interactions. Its stabilizing role on topological order has 
been investigated in terms of the topological entanglement entropy [54]. Here we are interested in 
how disorder enhances the reliability of a quantum memory. The underlying intuition is analogous to 
the reasoning used to explain why disorder suppresses zero-temperature electron transport in wires: 
the anyons, like electrons, get trapped in the many local minima present in a random potential. This 
phenomenon is commonly known as Anderson localization, named after the model introduced by 
Anderson to study the effect of randomness on quantum transport [5]. 

Perhaps the most well-studied paradigmatic model exhibiting Anderson localization is that of a 
single particle hopping on Z*^ with random on-site potential at each site. This corresponds to a 
tight-binding approximation of an electron moving in some random medium. The model Hamiltonian 
acting on states t/j G ^^(Z'^) takes the form Hji = —A + rjW, where A is the discrete Laplacian, 
{W'ijj){x) = W{x)tP{x) and where {W{x)}^i^^d are identically and independently distributed random 
variables with support e.g., on [—1,1]. The parameter ?? > controls the strength of the disorder in 
this setting. A realization of such a random Schrodinger operator is said to have exponentially 
localized eigenfunctions with a localization length if every eigenfunction satisfies ^ 



for some xq G Z*^ and some constant C > 1 depending on t/j. This is known to hold with high 
probability for arbitrary disorder strength in ID systems {d = 1). This fact can be established using 
transfer matrix methods and computing Lyapunov exponents [18, 19]. It is in stark contrast to the 
delocalized form of the eigenstates of the 'clean' Hamiltonian He = —A, which are Bloch waves. For 
d > 1, as well as for certain kinds of disorder, the situation is more complex: for example, there 
can be both localized and extended eigenstates separated by an energy threshold called a mobility 
edge. Theoretical manifestations of localization such as (1.1) or the so-called dynamical localization 
condition can be connected to experimentally measurable quantities, e.g., conductivity. Powerful 
analytic tools for establishing such conditions for more general families of random Hamiltonians have 
been developed (see e.g., [38, 52] for two recent reviews). This includes Hamiltonians with more general 
forms of disorder such as randomness in the hopping terms, i.e., off-diagonal disorder (see e.g., [20]). 

^It is worth mentioning that a typical reahzation of a random potential W{x) for an infinite lattice contains arbitrarily 
long intervals at which W{x) is arbitrarily close to a constant function. Such almost flat intervals would give rise to 
eigenfunctions which are delocalized at small scales. Note that the constant C in Eq. (1.1) could be arbitrarily large, so 
this bound does not prevent il^ix) from having oscillatory behavior for small \x — xo\. 
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A first step towards showing that Anderson localization increases the robustness of a topological 
qubit is to prove that the given physical system indeed has the desired localization properties. A 
difficulty here is that time-evolving topological memories generically involve coherent superpositions 
of states with varying number of anyons — a regime in which the standard Anderson localization 
theory cannot be directly applied. An approach pursued in [56, 51] is to restrict the dynamics to 
subspaces with a fixed particle number and neglect processes creating and destroying pairs of anyons. 
Within this approximation localization properties of the corresponding multi-particle states have been 
established both numerically [56] and analytically [51, 1, 16]. 

A second step is to connect localization properties of states to the robustness of the quantum 
memory. The first analysis along these lines starting from two-particle wavefunction localization was 
provided in [56, 51] for the toric code. By using the fact that logical errors correspond to topologically 
non-trivial trajectories of such pairs, some estimates on the failure probability starting from a two- 
particle configuration as well as certain geometrically arranged multi-particle initial configurations 
were obtained. However, these estimates neglect the full many-particle dynamics, and hence only 
provide a somewhat qualitative picture of the effect of localization. 

In order to assess the effect of disorder on the stability of topological memories more rigorously, 
we first need to choose a sufficiently simple model describing encoding, storage, and readout stages 
as well as some benchmark measure to compare the robustness of a memory in different regimes. 
Let us formulate this in more detail. We will describe the zero-temperature quantum information 
storage using a memory Hamiltonian Hq, a perturbation V, and an error correction operation $ec- 
The memory Hamiltonian Hq = — Sj represents an ideal topological quantum memory such as 
the toric code model. It is a sum of local pairwise commuting stabilizers Sj where each stabilizer 
has eigenvalues ±1. Ground states of Hq are +1 eigenvectors of any stabilizer Sj while excitations 
correspond to flipped eigenvalues, Sj = —1. The encoding stage amounts to initializing the system in 
some ground state \g) of Hq, see Section 2.2. The ground state of Hq must be degenerate to permit 
encoding of one or several qubits. The storage is modeled by a unitary evolution under a Hamiltonian 
H = Hq + V for some time t which results in a final state \g{t)) = e*(^"+^)* \g). Here we assume 
that the perturbation V is either unknown or known only partially, so the unitary evolution cannot 
be undone even if we have a full control of the system. However, if the perturbation V is sufficiently 
weak and the evolution time t is not too long, the final state \g{t)) can be regarded as a slightly 
corrupted version of the initial state 1^) with a small density of errors which can be corrected at 
the readout stage. The readout begins by measuring a syndrome, i.e., measuring the eigenvalue ±1 
for every stabilizer Sj. The syndrome information s is fed into an error correction algorithm that 
determines a unitary correcting operator C(s) returning the system back to some ground state of 
Hq. The error-corrected final state is therefore ^eciloit)) {g{t)\) , where <I>ec is a linear map describing 
how the error correction acts on states. More explicitly, ^ec(p) = '^s^(^)QsPQsC{s)\ where Qg is 
the projector onto a subspace with a fixed syndrome s, see Section 2.3 for details. We will measure 
robustness of a memory by its storage fidelity — the overlap between the initial encoded state and the 
final error-corrected state: 

= (5l<J>ec(e'(^°+'')*|5)(5|e-'(^°+^)*)|5>, (1-2) 

and the worst-case fidelity F{t) = min|g^ i^^^ (t) minimized over all encoded states. By definition, 
< F{t) < 1 with F{0) = 1. A perfect recovery of encoded information corresponds to F{t) = 1. 
Another measure we will use is the storage time Tgtorage = 7storage(-Pb)- This is the time it takes 
for the storage fidelity F{t) to drop below a given threshold value Fq, say, Fq = 0.99 (we note that 
although F{t) tends to decrease with increasing t, its behavior may be highly non-monotone, see e.g. 
Fig. 3(b) in Section 7). We note that error correction is necessary here to give a fair assessment of 
the information recoverable from the memory. But its role is different from the case of active error- 
correction: it is only applied once in this process when the information is retrieved after time t. In 
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particular, it plays no direct role in preserving coherence in the time interval [0,t]. 

Quenched disorder will be modeled by introducing some randomness into the perturbation V. 
Clearly, physically realistic models of disorder depend on the system under consideration. Here we 
study the Majorana chain model proposed by Kitaev [40]. It describes a quantum wire put in contact 
with a bulk superconductor which enables creation and destruction of electron pairs in the wire through 
tunneling of Cooper pairs (see Section 2.1 for a formal definition of the model and a discussion of its 
properties). Since this is a system of electrons, the simplest choice of disorder is a random site- 
dependent chemical potential. The corresponding perturbation is F = —^f=i fJ'jCijCLj, where (oj) 
are electron creation (annihilation) operators, N is the system size, and fij is the chemical potential 
at the site j. Furthermore, we assume that fij = fj, + rjxj, where Xj G [— 1,1] are independent 
identically distributed random variables drawn from the uniform distribution, while rj > controls the 
disorder strength. The constant /U represents the homogeneous part of the perturbation. More general 
perturbations that might be relevant for the Majorana chain model are described in Section 2.1. We 
note that our description of quenched disorder differs from the one of [56, 51] where the randomness 
was introduced into the memory Hamiltonian Hq by allowing site-dependent random coefficients in 
front of the stabilizers Si. The reason is that the latter choice of randomness would be very artificial 
and hard to motivate in the case of Majorana chains, see Section 2. 

In the presence of disorder the storage fidelity F{t) is a random variable depending on the disorder 
realization {xj}. We will mostly be interested in the expectation IE[i^(i)] over disorder realizations. 
Comparing the behavior of F{t) for the clean case (?] = 0) and E[F(t)] for the disordered case {rj > 0) 
gives a precise way of studying the effect of disorder on the robustness of quantum memories. Here 
we are particularly interested in the asymptotic scaling of the expected storage time E[rstorage] in the 
limit of large system size. 

Main results 

Our results are three-fold. First, we formulate a dynamical localization condition for the one-particle 
Hamiltonian which is sufficient to achieve an exponential scaling of the expected storage time, that is, 

E[rstorage]=e^(^), (1.3) 

in the limit of large system size N, see Theorem 2 in Section 2.4 for details. We conjecture that our 
dynamical localization condition is satisfied in the limit of weak perturbations and strong disorder, that 
is, <C <C 1. In other words, the homogeneous part of the perturbation must be small compared with 
the random part. We prove that the required localization condition is satisfied whenever the entries 
of the orthogonal matrix describing the time evolution of the Majorana modes decay exponentially 
away from the diagonal with an A^-independent localization length ^i, see Lemma 3 in Section 6.3. 
The localization length may diverge in the limit 77 — )• 0, as is the case in the standard ID Anderson 
model, but it must be upper bounded as ^1 = 0{r]~^) for some sufficiently small constant 7 > 0. 

Second, we give supporting evidence that the desired scaling of the localization length can be 
achieved by computing Lyapunov exponents of the one-particle eigenfunctions, see Section 6.2. This 
suggests a scaling ~ log(r/~^) in the limit 77 — )• when the ratio fi/r] is kept constant. We note 
that the logarithmic divergence of the localization length at the band center is a common feature of 
systems with disorder in the hopping amplitudes (so called off-diagonal disorder), see e.g. [22, 18]. 

Third, we compute the storage time numerically using the Monte Carlo technique developed by 
Terhal and DiVincenzo [53] for simulating quantum dynamics and measurements for non-interacting 
fermions (see Section 7). The running time of our algorithm grows as where 6 is the precision 

up to which one needs to estimate the storage fidelity. It allows us to compute the storage time for 
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chains with a few hundred sites (up to = 256) in the regime of strong perturbations^, that is, /x ~ 1 
and 1] = (clean case), and r/ ~ /x ~ 1 (disordered case). The simulation shows that in the absence of 
disorder the storage time grows as a logarithm of the system size: 

Tstorage~0(logA^). (1.4) 

This scaling has been recently predicted by Kay [36] based on mean-field arguments, see Section 7 for 
details. In the presence of disorder we observe an approximately linear scaling E[Tgtorage] ~ This 
confirms the expected enhancement of the storage time, although the enhancement is much weaker 
than our theory predicts, see Eq. (1.3). This discrepancy could be accounted for by the fact that 
the system size is comparable with the localization length in the simulated regime, whereas 
Eq. (1.3) is expected to hold only when N ^ ^g. It could also point to an interesting possibility that 
a crossover from a polynomial to an exponential scaling of the storage time occurs as one interpolates 
between strong (// ~ ry ~ 1) and weak (^ <C <C 1) perturbations. Finally, we give examples where 
an artificially engineered deterministic disorder potential leads to improved storage times compared 
to random disorder. 

For a more formal statement of our results see Section 2.4. 

Discussion 

The Majorana chain model [40] has a number of properties which make it amenable to the study of 
localization. A major difference to e.g., the 2D toric code is that the perturbed memory Hamilto- 
nian Hq + is an exactly solvable model, assuming that the perturbation V is quadratic in fermionic 
operators. This means that the full many-particle dynamics can be connected to a single-particle prob- 
lem. The exact solvability of the model provides powerful analytical tools for analyzing the effect of a 
perturbation, see e.g. Theorem 3 in Section 2.4 and makes possible efficient numerical computation of 
the storage fidelity. Finally, Anderson localization in the ID geometry requires only arbitrarily weak 
disorder which is likely to be present in the system due to natural imperfections (see Section 6.2). 

Practical proposals for physically realizing the Majorana fermion chain model have recently re- 
ceived a lot of attention. Fu and Kane [24] showed that a ID wire with Majorana edge modes can be 
realized in a superconductor-topological insulator-superconductor junction. This exploits the prox- 
imity effect between an ordinary (s-wave) superconductor and the surface of a strong topological 
insulator. Several groups proposed to realize and compute with Majorana fermions in ID semicon- 
ducting wires or networks deposited on an s-wave superconductor [50, 44, 2]. More recently, Jiang et 
al. [35] suggested a realization of Majorana fermions using optically trapped fermionic atoms. Methods 
for accessing the encoded information were studied in [34, 6, 33, 29]. 

One may wonder whether the paradigmatic toy model of a disordered memory studied here bears 
any relevance to these practical proposals. Especially with realizations based on superconductors, 
a major difference to our idealized model is that most physical operations directly couple to the 
topological charge, i.e., the quasi-particles defined by the system Hamiltonian. This is fortunate 
if we consider a regime of small system size (compared to the inverse perturbation strength) and 
short storage times: here the information may be thought of as evolving under a Hamiltonian without 
hopping terms, rendering error correction unnecessary (cf. Corollaries 1,2 in Section 2.4 and Fig. 3(a) in 
Section 7). However, outside this regime, error-correction is essential to preserve information for longer 
times. Here the fact that operations couple to quasi-particle excitations is less welcome: it implies 
that encoding and measuring syndrome information corresponding to the 'clean' system as envisioned 
here may be difficult or impossible. Ultimately, this depends on the system under consideration. Here 

^As we discuss below, the storage fidelity is close to 1 whenever e < simply because the perturbed ground 

state has a large overlap with the unperturbed one. To explore the asymptotic scaling of the storage time for weak 
perturbations, say, e ~ 10"'^, one would need to simulate chains with at least ~ 10'' sites. 
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we take the viewpoint that some form of error-correction will most likely be required in any realization 
of a stable quantum memory. This justifies the consideration of idealized error-correction procedures, 
but also poses the challenge of identifying systems where such operations may be practically realizable. 

Outline 

In Section 2, we introduce Kitaev's model [40] of quantum wires with unpaired Majorana modes, 
describe the error-correction procedure and state our main results. In Sections 3,4 we prove two 
technical theorems needed for our analysis of the storage fidelity. Theorem 1 provides a tight bound 
on the strength of a perturbation capable of destroying the unpaired Majorana modes. Theorem 3 
is a stronger version of the quasi- adiabatic continuation method [45, 31, 30] in which the effective 
Hamiltonian that governs the quasi-adiabatic evolution has a strictly exponential decay of interactions. 
These theorems are used in Section 5 to show that the storage time scales exponentially with the system 
size in the Anderson localization regime (Theorem 2). In Section 6, we discuss localization properties 
of the one-particle Hamiltonian and compute Lyapunov exponents. Our simulation algorithm and 
numerical results are discussed in Section 7. 



2 Quantum wires with unpaired Majorana fermions 
2.1 Definition of tiie model 

Following [40], we consider a system of spinless electrons that live on a chain with N sites. The 
electrons can hop between adjacent sites. In order for unpaired Majorana modes to emerge, the chain 
must be put in contact with a bulk superconductor. In this setup the number of electrons in the chain 
is only conserved modulo two due to the tunneling of Cooper pairs between the bulk superconductor 
and the chain. This model can be described by the following Hamiltonian proposed in [40]: 

N N-1 

H = -^j{a}-aj - 1/2) + ^ -w{a)-aj+i + a^j^^aj) + Aa^aj+i + A*aJ^^aj, (2.1) 
i=i i=i 

where is the chemical potential at a site j, is the hopping amplitude, and A is the superconducting 
pairing induced by the proximity effect. Note that A is typically complex, A = |A|e*^. We will refer 
to the Hamiltonian (2.1) as a Majorana chain model. It will be convenient to rewrite H in terms of 
2N Majorana fermionic modes ci, . . . , C2n by 

aj = — ^ — (c2j-i + ic2j) and a] = {c2j~i - ic2j) , j = l,...,N. (2.2) 

The Majorana operators are self-adjoint, Cp = Cp, and obey commutation rules 

CpCq + CqCp — 25p ql , Cp — / , (^•3) 

where I is the identity. We can now rewrite the Hamiltonian Eq. (2.1) as 

H = Ho + V , (2.4) 

where 

iJo = Y ^ C2jC2j+i , J = A + w (2.5) 
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will be chosen as our 'clean' memory Hamiltonian and 



V 



lAI 



Af-l 



W 



lC2j+2 



TV 



-lC2j 



(2.6) 



We will be mostly interested in the regime |A| ~ w ^ so that Hq corresponds to the exact 
equalities w = |A|, /Uj = 0, while V represents small deviations from this idealized setting. It will 
always be assumed that the parameters A and w are known, while the chemical potential plays 
the role of an unknown perturbation. The encoding and decoding steps should not depend on the 
realization of jjLj. In our simulations the chemical potential will be chosen as 



X, 



e [-1,1], 



(2.7) 



where xi, . . . ,xn is a realization of the quenched disorder and rj > controls the disorder strength. 
We will consider three scenarios (in all cases we assume ^ > 0): 



No disorder 


r? = 0. 


Random disorder 


rj > and {xj} are i.i.d. random variables with uniform 


probability density. 


Pseudo-random disorder 


rj > and the sequence {xj} is optimized to minimize the 
localization length. 



In the first case (no disorder) the system is translation-invariant and excitations of H can be 
described by Bloch waves'^. This clean case will provide us with baseline estimates of the storage 
time needed to assess the effect of adding randomness. The second case (random disorder) attempts 
to describe a random external electric potential, for example, the one created by impurities. In this 
regime coherent propagation of quasiparticles is suppressed due to the Anderson localization and a 
more favorable scaling of the storage time is to be expected. Once it is established that Anderson 
localization is beneficial for the memory stability, a natural question is whether it can be brought 
about on purpose by applying a suitably engineered external potential. This is the main motivation 
for considering the third case (pseudo-random disorder). Here we assume that the constant part of 
the potential // represents an unknown perturbation, while the sequence {xj} represents a controlled 
external potential introduced on purpose in order to localize wavefunctions of the quasiparticles. For 
simplicity we set w = |A| in the simulations, while our analytical results summarized in Section 2.4 
hold for any w and A. 



2.2 Encoding a qubit into the ground state 

The Hamiltonian Hq is a sum of pairwise-commuting terms formed by disjoint pairs of modes C2jC2j+i 
which can be diagonalized simultaneously. Every pair of modes C2j, C2j+i can be formally combined 
to a complex fermionic mode with an excitation energy J. Ground states of Hq are —1 eigenvectors 
of every term ic2jC2j+i which can thus be regarded as 'stabilizers' of a quantum code. The boundary 
Majorana modes ci and C2n play ^ special role since they do not appear in Hq and commute with 

^The analogue of Bloch waves for a chain with boundary are superpositions of right-moving and left-moving waves. 
The relative amplitudes of these waves and the admissible values of the momentum can be easily computed using the 
scattering method. 
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every term in Hq. Combined together, ci and C2 form a complex boundary mode whose occupation 
number operator is 

hb = b^b, b = ^{ci + ic2N) ■ (2.8) 

This boundary mode has zero excitation energy since [6, -f^o] = 0. Let Ig^j) be the ground state of Hq 
in which the boundary mode has occupation number hf, = a, o" = 0, 1. 

Although the ground state of Hq has two-fold degeneracy, it cannot be used directly to encode 
a qubit. The reason is that {go) and l^i) belong to different superselection sectors defined by the 
fermionic parity operator 

N 

p = (_i)Ef=i4% = l[{-i)c2j-iC2j . (2.9) 

i=i 

Rewriting P as a product of C2jC2j+i and ciC2n one can easily check that P \ga) = {—iy\ga)- Since any 
physical state has either an even or odd number of fermions, we conclude that coherent superpositions 
of \go) and l^ii) are unphysical. We will avoid this problem by adding a reference system R that 
contains one or several fermionic modes with some special states |0/j) and having even and odd 
fermionic parity. We can then encode a one-qubit state a|0) + /3|1) into a state 

|5) =a|5o)»|Ofl) (2.10) 

in which the total number of fermions is even. For simplicity we will assume that the reference 
system has trivial dynamics, although our results can be generalized in a straightforward way to the 
case when the reference system is another Majorana chain with the Hamiltonian Eq. (2.1), that is, 
|0i?) = \go) and \1r) = \g\). In this case we can choose logical Pauli operators on the encoded qubit 
as Z = [—i)ciC2N and X = (— i)cici, where ci is the first Majorana mode on the chain representing 
the reference system. 

2.3 Error correction 

Let us now discuss how one can perform the decoding, that is, how to retrieve the original encoded 
state l^f), see Eq. (2.10), from the time-evolved state \g{t)) = e^^^\g). Note that we cannot simply 
reverse the time evolution since some terms in H are unknown. 

The important insight made in [40] is that the ground subspace spanned by l^o) and l^i) can be re- 
garded as a quantum error correcting code protecting encoded information against local perturbations 
that involve even fermionic operators, i.e., operators commuting with P. Let us first characterize the 
set of errors that can be corrected by this code. Define stabilizers 

Sj = {-i)c2jC2j+i , j = 1, . . . , - 1 , (2.11) 
such that Hq = — J2f=i^ ^^^"^ elementary errors 

Ej = = {-i)c2j-iC2j , j = 1, . . . , A (2.12) 

that are analogous to single-qubit errors in the standard error correction theory. The unitary evolution 
^iHt ^g^^ represented as a linear combination of errors of a form E = Cj^ ■ ■ ■ Cj^^ , where each error E 
can be uniquely written as a product of stabilizers and elementary errors (up to an overall phase). Any 
such error either commutes or anti-commutes with any stabilizer Sj. We will say that E is an error 
of weight V if it involves exactly v elementary errors. One can easily check that an error E commutes 
with all stabilizers only if f = or = A. In the first case £^ is a product of stabilizers and thus it 
acts trivially on any ground state \g). In the second case E \s a, product of stabilizers and the parity 
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operator P, see Eq. (2.9), which can be regarded as the logical-Z operator since P\ga) = (— 1)'^|5'ct)- 
Note that P is equivalent to {—i)ciC2N modulo stabilizers. 

Error correction begins with the syndrome measurement. It involves a non-destructive eigenvalue 
measurement for every stabilizer Sj which yields an eigenvalue (— 1)*^ , Sj G {0, 1}. Let s G {0, 
be the syndrome, that is, the list of all measurement outcomes. The syndrome measurement maps the 
state \g{t)) to Qs\g{t))-, where 

N-l 

n 2 (^+^"^^''^^) ^^-^^^ 

is the projector onto the subspace with syndrome s. (Note that the stabilizers Sj depend on the 
superconducting phase 6, namely, Sj = ajcj+i + a^jj^idj — e^^Oj-aj+i — e~*^aj_(_]^aj. Hence the syndrome 
measurement can be realized only if the superconducting phase 9 is known.) 

We can now return the system back to the ground state using the syndrome information s by 
applying a correction operator C{s). It will be chosen as a product of elementary errors Ej which is 
(a) consistent with the observed syndrome, and (b) uses as few elementary errors as possible. The 
first condition demands that 

C{s)Sj = {-iy^SjC{s) (2.14) 

for all j = 1, . . . , — 1. It guarantees that the corrected state C{s)Qs\g{t)) is indeed in the ground 
subspace of Hq since it is a +1 eigenvector of every stabilizer Sj. The second condition captures the 
intuition that the Hamiltonian H involves only errors of small weight {v = 1 and v = 2) and thus typical 
errors generated by the unitary evolution e*^* have weight much smaller than N (which may or may not 
be true depending on the evolution time t). If E' and E" are products of elementary errors consistent 
with the syndrome s then the product E'E" commutes with all stabilizers and thus either E' = E" or 
E' = E"P = E" Y\f=i Ej (modulo stabilizers). It follows that Eq. (2.14) has only two solutions which 
use mutually complementary subsets of elementary errors. The smallest weight solution C{s) thus has 
weight at most N/2. We note that the error correction can only fail by introducing a logical-Z error 
since it uses only even fermionic operators. In other words, C{s)Qs\g{t)) ~ (og/ + bsP)\g) for any 
observed syndrome s with some amplitudes as,bs G C. Also note that in contrast to stabilizers, the 

elementary errors Ej = (— 1)°j"j and the correction operators C{s) do not depend on the parameters 
of the Hamiltonian. 

The procedure described above is analogous to the minimum weight decoding for the ID repetition 
code, see for instance [21]. If a syndrome s is viewed as a collection of excitations, a candidate error 
E that might have caused s can be viewed as a pairing between the excitations. Every excitation in s 
is paired with some other excitation or with the boundary by a chain of elementary errors Ej. Since 
in the ID geometry any set of excitations can be paired in only two different ways, any syndrome s 
has exactly two candidate errors E',E" consistent with s (modulo stabilizers) such that E'E" ~ P 
(modulo stabilizers). The correction operator C(s) is chosen as the shortest error chain. 

One can formally describe how the error correction acts on states by a TPCP map 

^ec(/>) = C{s)QspQsC{s)\ (2.15) 

s 

where the sum is over all 2^~^ syndromes s. This gives the following expression for the storage fidelity 
with an initial state \g) and a storage time t: 

F|,)(t) = {g\^Mt)){9{t)\)\g) = \{9\C{s)Qse'%)f . (2.16) 

s 
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Since all operators H,Qs, and C(s) preserve fermionic parity, the reference system plays no role in 
the expression for i^g^(t). In particular, one can formally compute F\g){t) using an 'unphysical' initial 
state \g) = a\go) + P\gi) instead of the one defined in Eq. (2.10). 

It is worth pointing out that, as discussed in the introduction, error correction is not the only 
possible strategy to protect the encoded qubit from decoherence. In the special case of the memory 
Hamiltonian and the perturbation described by Eqs. (2.5,2.6) the refocusing technique can be used 
instead. For example, applying a refocusing operator C = C1C3 . . . C2Ar-i exactly in the middle of the 
time evolution would reverse the sign of the Hamiltonian making the overall evolution trivial: 

(j\^iHt/2(j^iHt/2 ^ ^-iHt/2^iHt/2 ^ j _ 

The advantage of the error correction approach is that it applies to perturbations much more general 
than the one defined in Eqs. (2.6), see below. 

2.4 Summary of main results 

It will be convenient to state our results for perturbations more general than the one defined in 
Eq. (2.6). Specifically, we will assume that = (V^) X]j=i^ C2jC2j+i is the Hamiltonian defined in 
Eq. (2.5) with J = 1 and 

. 2N 

^ = i E ^f.^^P^^- (2.17) 

p,q=l 

Here V is some real anti-symmetric matrix of size 2N. Let us say that a perturbation has strength e 
and range r iff 

||F|| < e and Vp^g = unless \p - q\ < r . (2.18) 

For example, the perturbation V describing the Majorana chain, see Eq. (2.6), has strength e = 
maxj l/iji + \w — |A|| and range r = 3. Note that any physical perturbation V has strength 0(1) even 
though the norm of ||y|| is typically an extensive quantity, that is, ||V^|| ~ A^. 

Recall that the unperturbed Hamiltonian Hq has a two-fold degenerate ground space separated 
from excited states by an energy gap AE = J = 1. The degeneracy stems from two unpaired Majorana 
modes localized on the left and the right boundaries of the chain. Our first result provides a tight 
bound on the perturbation strength capable of destroying these zero-energy boundary modes and 
closing the gap above the ground state. This result can be regarded as a stronger version of the gap 
stability proved for more general Hamiltonians with topological order in [9, 8]. Unfortunately, the 
techniques used in the present paper only apply to non-interacting fermions. 

Theorem 1 (Gap stability under perturbations). Consider any perturbation V with strength 
e < 1 and range r. Let < E'l < E2 be the three lowest eigenvalues of Hq + V . Then 

eI-EI<[-^)" and 4-4>i_,_(^Ji_y'- (2.19) 
for all large enough N . 

Hence any short-range perturbation with strength e < 1 preserves the two-fold ground state de- 
generacy of Hq up to exponentially small corrections. The degenerate ground state is separated from 
excited states by an energy gap roughly 1 — e. To see that a perturbation with strength e = 1 is capable 
of closing the gap E^ — e\, we can choose V = {ie/2) Ylf=i C2j-iC2j = ^Ylf=iiO']^j ~ l/^)- In this 
homogeneous case the gap of .^0 + ^ was calculated in [40] , namely, E2 — e'I = 1 — e ior all < e < I. 
In that sense the bound of Theorem 1 is tight. Let us emphasize that the energy splitting E^ — Eq for 
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the ground state is bounded by a pure exponential function of A^. In contrast, more general bounds 
on the energy splitting proved in [9, 8] decay slower than exponentially. The proof of Theorem 1 also 
provides an explicit construction of the unpaired Majorana modes corresponding to the perturbed 
Hamiltonian Hq + V. These modes are linear combinations of the Majorana operators cj which are 
exponentially localized on the left and on the right boundaries of the chain 
Our second result is a lower bound on the worst-case storage fidelity 

F{t) = rnin (5|^ec(e^(^»+^)*b) (<?|e-*(^°+^)*)b) , (2.20) 
b) 

or more precisely, its expectation E[F(t)] over disorder realizations. We give a dynamical localization 
condition for the one-particle problem which, when satisfied, implies that E[F(f)] remains close to 1 
for times t exponentially large in the system size. 

Define an orthogonal matrix R{t) G S0{2N) describing evolution of the Majorana operators Cp in 
the Heisenberg picture, 

2N 

cpit) = e^(^"+^)*Cpe-(^«+^)* = ^ i?p,,(t)c, . (2.21) 

q=l 

In the case when V includes random disorder, the theory of Anderson localization roughly predicts 
that the time-evolved operator Cp{t) is well-localized near the original mode Cp for all t G M with high 
probability over the disorder distribution. This is equivalent to the so called dynamical localization 
condition: the matrix of expectation values E[|i?p^g(t)|] decays rapidly away from the main diagonal 
(see e.g., [52]). Our bound on the storage fidelity requires dynamical localization only at large scales, 
that is, |p — g| ~ A^. Unfortunately, we also need to impose constraints on the decay of certain multi- 
point correlation functions composed of matrix elements of i? = R{t). More precisely, suppose p,q 
are m-tuples of integers in the interval [1, 2A^] such that pi < p2 < ■ ■ ■ < Pm and qi < q2 < . . . < qm- 
Define a distance between p and q as 

m 

a=l 

We denote by q\ the mx m submatrix of R obtained by retaining only the rows indexed by p and 
the columns indexed by q, i.e., R\p, q]a,b = RpaSb all 1 < a, 6 < m. 

Definition 1 (Multi-point dynamical localization). The unitary evolution operator ^^{^o+V)t 
exhibits multi-point dynamical localization iff there exist constants C, ^ > such that for all sufficiently 
large N , for all m > 1, and for all ordered m-tuples p, q such that \p — q\i > one has 

E[\det R[p,q]{t)\] < C'"e-^/«. (2.22) 

Here the expectation value is taken over disorder realizations. 

We will refer to the constant ^ as the (dynamical) localization length. In general, ^ does not coincide 
with the localization length that controls decay of eigenfunctions (the spectral localization). The 
relationship between ^ and is discussed in more detail in Section 6. For m = 1 one has det R[p, q] = 
Rp,q and the multi-point localization condition reduces to the usual (two-point) dynamical localization, 
although we only need it for large scales, |p — (?| ~ N . In Section 6.3 we discuss how to replace the 
expectation value of the determinant in Eq. (2.22) by more standard multi-point correlation functions. 
Let us emphasize that in the present paper we do not prove Eq. (2.22). Instead, our goal is to establish 
a link between dynamical localization and the storage fidelity. It is provided by the following theorem. 
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Theorem 2 (Storage fidelity). Consider any random ensemble of perturbations V with strength 
e < 1/4 and range r. Consider any time t such that the unitary evolution g^i^o+V)t gf^^yg ^/^g fnulti- 
point dynamical localization condition with some localization length ^ and constant C . Then the 
expected storage fidelity can be bounded as 

E[F(t)] > 1 - t2 . e-v - 0(l) . e-*^^ (2.23) 

for all sufficiently large N , where 



^ = ^' ^' ^) = max(g,16r,280) " "^'^ " ' 
and a, /3 > are some constants that depend on details of the proof. 

We conclude that for a fixed fidelity threshold Fq < 1, the storage time rstoragc(^o) grows expo- 
nentially with N if the ensemble of Hamiltonians satisfies multi-point dynamical localization with 
time-independent constants C, such that u > {). Computing the parameters i^, C and checking 
whether v > Q obviously requires more detailed information about the perturbation V and the disor- 
der distribution. In the special case of the Majorana chain model we have e < ^ + rj, where yu and rj 
control the strength of the homogeneous and the random parts of the perturbation, see Eq. (2.7). In 
Section 6.2, we provide some evidence indicating that the regime z/ > can be achieved for /i ~ r/ ^ 1 
by computing the Lyapunov exponent i of the one-particle eigenfunctions. With the technique de- 
veloped by Eggarter et al [22] for analyzing tight-binding chains with an off-diagonal disorder we 
show that ~ log (1/r/) in the limit 77 — when the ratio ///r/ is kept constant. Using the rough 
estimate ~ ~ and assuming C = 0(1) in this limit, we would get 

z.(e, C, r) « i - (3Ce^e = - 0(1) • r?i^ > 

4 log(l/??) 

for small enough rj. Another scenario where the condition > could be satisfied is when some 
randomness is present in the coefficients of the unperturbed Hamiltonian Hq so that the disorder 
strength and the perturbation strength are independent parameters. In this scenario we expect bounds 
of the form ^, C = 0(1) in the limit e — >• 0, so that u > for sufficiently small e. Our proof of Theorem 2 
can be easily extended to this setting. 

Our last theorem is a technical tool needed to prove Theorem 2, although it might be interesting 
on its own right. It provides a stronger version of the exact quasi-adiabatic continuation [45, 31, 30], 
that is, a unitary operator mapping the ground subspace of Hq to the one of Hq + V . 

Theorem 3 (Quasi-adiabatic continuation). Consider any perturbation V with strength e < 1/2 
and range r. Let \g(j) and be the ground states of Hq and Hq + V in the sector with fermionic 
parity a G {0,1}. Then \ga) = U ■ \ga) , where U is a unitary operator describing evolution under a 
time- dependent quadratic Hamiltonian D[u), u G [0,1], with strength roughly e. More precisely, 

1 



C/ = r-exp i duD{u) , (2.24) 
L ^0 

where 

. 2N 

D{u)='- DpA^)^P^'i (2-25) 

p,q=l 

for some real anti- symmetric matrix D[u) such that 
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for all u G [0, 1] and all N . Furthermore, i/e < 1/4, the Hamiltonian D{u) is exponentially decaying: 

\Dp,qiu)\ < cer^ (2/3)^ for allp,q . (2.27) 

The constant coefficient c = 0(1) depends on details of the proof. 

Note that this theorem requires a stronger bound on e although we beheve that the exponential 
decay of Dp^g{u) holds for all e < 1. Our proof yields a rather large constant coefficient c = 2 x 10^ 
since we have not tried to optimize it. Let us also note that a weaker version of Theorem 3 in which 
\Dp^g{u)\ has stretched exponential decay can be easily derived using the techniques of [45, 31, 30]. 
Unfortunately, any bound in Eq. (2.27) decaying slower than exponentially is not sufficient for our 
purposes since it results in an error distribution in which uncorrectable high-weight errors are not 
sufficiently suppressed. 

A simple corollary of Theorem 3 is that the overlap |(ffo-|(7o-)| between the unperturbed and per- 
turbed ground states is close to 1 as long as N <^ e~^. 

Corollary 1 (Perturbed versus unperturbed ground states). Consider any perturbation V 
with strength e < 1/8. Let \gu) and Ig^j) be the ground states of Hq and Hq + V with fermionic 
parity a £ {0, 1}. Then 

K5<xb<x)l' > l-4iVe2 (2.28) 

for all (T = 0, 1 and for all N. 

This also implies that error correction is not really necessary in the regime <C e^^, since the 
initial encoded state \g) is very close to the ground state of Hq + V. As a consequence, the unitary 
evolution e^iH'^+V)^ has no effect on \g) except for an exponentially small dephasing between Igo) and 
\gi). For simplicity, below we consider the maximally entangled encoded state, \g) = (j^o) \0r) + 
\gi) ^ |li?))/\/2, see Section 2.2. The quantity F\g){t) is the entanglement fidelity [32] of the channel 
composed of storage, time-evolution and error-correction. 

Corollary 2 (Cosine law). Consider any perturbation V with strength e < 1/8 and let 6 = — Eq 

be the energy splitting of the ground state, see Theorem 1. Then 

\F\g){t) - cos^ {6t/2)\ < 8e\/iV . (2.29) 
This bound applies even if no error correction is performed, that is, if <I>ec is the identity map. 
We conclude that the interesting region of parameters in which error correction really matters is 

Finally, let us briefly summarize the results of numerical simulations performed for the Major ana 
chain model with w = |A| and a random site-dependent chemical potential fij = fi + rjxj, where 
Xj G [—1, 1] are i.i.d. random variables, see Section 7 for details. The corresponding perturbation 
strength is e < /j, + r]. We computed the storage fidelity and the storage time in the following three 
regimes: 

Weak perturbations (e <^ 1): In this regime we can only perform the simulation for A^ <^ l/^^- 
Here our numerical results agree with Corollary 2: the storage time is solely determined by the 
dephasing of the ground states and grows exponentially with the system size, Tgtorage ~ 1/5 = 
exp(r2(A^)), where 5 is the ground state energy splitting, see Theorem 1. Disorder does not play 
any role here. 

Strong perturbations (e ~ 1), no disorder: Here we observe a logarithmic scaling of the storage 
time, Tstorage ~ log(A^), for 4 < A^ < 256. This behavior can be understood by mapping the 
time evolution of the Majorana chain model to a quantum quench in the transverse-field Ising 
chain and using known results on the decay of magnetization in the latter model, see Section 7. 
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Strong perturbations with disorder: We observe an approximately linear scaling of the storage 
time, IE[Tstorage] ~ for 4 < < 128. The localization length is comparable with the system 
size for this range of A^. The storage times are clearly enhanced compared to the clean case. At 
the same time, we expect that the exponential scaling of the storage time cannot be achieved for 
strong perturbations, e ~ 1, because < in Theorem 2. The storage fidelity exhibits strong 
fluctuations as a function of disorder; no self- averaging behavior is observed. 

Surprisingly, we also found that the random i.i.d. disorder potential {xj} is not the optimal one as 
far as the storage time is concerned. The best performance was numerically observed for a disorder 
potential generated by iterations of the logistic map, namely, xj = 1 — 2yj, where y^+i = ayj{l — yj) 
for 1 < J < A'^ — 1. Such disorder potentials are uniquely determined by the initial condition yi £ [0, 1] 
and the parameter a. Note that yj £ [0, 1] for all j as long as < a < 4. Sequences generated by 
iterations of the logistic map are known to exhibit chaotic behavior for 3.57 ^ a < 4 and almost all 
initial conditions. For example, choosing a = 3.9914, yi = 0.2845 we were able to boost the storage 
time from E[Tstoragc] ~ 31 to Tstoragc ~ 176 for A^ = 64. 



3 Unpaired Majorana modes are stable under perturbations 

The proof of Theorem 1 is presented in Section 3.2. We begin with a summary of the relevant 
connections between the one-particle and the many-particle problem in Section 3.1. 



3.1 Normal form of quadratic fermionic Hamiltonians 

We will be using the following well-known fact throughout the paper. 

Fact 1 (Williamson eigenvalues). Let A be any real anti- symmetric matrix of size 2N . Then there 
exist non-negative real numbers < Ai < . . . < Aat called Williamson eigenvalues of A and complex 
vectors ipi, . . . , ipN E C^^ such that 

Aipj = iXjipj and Aipj = —iXjipj 

where tpj is the complex conjugate ofipj. In addition, {ipj\ipk) = ^j,k o.nd {ipj\ipk) = 0. 
Consider a Hamiltonian 

2N 



p,q=l 



where H is an anti-symmetric real matrix. Then the Williamson eigenvalues 0<Ai<...<Aivofff 
are the single-particle excitation energies of H. To see this explicitly, let {V'jlj be the corresponding 
complex (normalized) eigenvectors, and let {\p)} be the computational basis of C^^. Fact 1 implies 
that the operators {aj}jLi defined by 



27V 



^ p=l 



obey the fermionic CCRs, that is, we can regard them as annihilation operators of some fermionic 
modes. A simple calculation shows that 

[H , aj] = — Ajdj , 
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hence the {aj} are canonical (complex) modes of H. In particular, if Eq is the ground state energy of 
H, then 

N 

H = El + Y,^ja]aj . (3.1) 
j=i 

Alternatively, we can write the Hamiltonian in terms of Majorana modes {cp}^^^ as 

. N 

6 = -^i^-i^aj where C2j-i = aj + a] and C2j = (-«)(aj - a^) • (3.2) 

The ground state energy Eq = Si=i '^^^ found, e.g., by comparing (3.1) with (3.2). 

3.2 Stability analysis based on the Brillouin-Wigner method 

Define anti-symmetric real matrices Hq and V such that 

. 2N . 2N 

i^o = ^ ^ {Ho)p,gCpCg and = - ^ Vp^gCpCg. (3.3) 
p,q=l p,q=i 

Let < Ai(n) < . . . < XNiu) be the Williamson eigenvalues of Hq + uV, where u G [0, 1]. Equivalently, 
Xj{u) are single-particle excitations energies of the Hamiltonian Hq + uV. Our choice of Hq implies 
that 

Ai(0) = and Xj{0) = 1 for all j = 2, . . . , . (3.4) 
Since \\V\\ < e, Weyl's inequality [55] yields 

Ai(u) < e and Xj{u) > 1 - e for all j = 2, . . . , TV . (3.5) 

In the rest of the proof we focus on the (normalized) eigenvector ipi = ipi (u) such that 

(Hq + uV)^i = iXiiu)i;u (V'i|Vi) = l, u£[0,l]. (3.6) 

Define ^ 

u* = max u subject to Xi{v) < for all v £ [0,u] . (3.7) 

uG[0,l] 2 

Continuity of Xi{u) implies that either u* = 1 or u* £ (0, 1) in which case Xi{u*) = (1 — e)/2. We 
will show that the second case actually leads to a contradiction for sufficiently large N (this trick is 
borrowed from [8]). From now on we will only consider u G [0,ti*], so with Eq. (3.5), we have 



Xi{u) < -(1 - e) and Xj{u) > 1 - e for ah j = 2, . . . , iV . (3.8) 



We will now explicitly compute ipi using Brillouin-Wigner perturbation theory. Let p = 

1, . . . , 2A^, be the standard orthonormal basis of in which Hq and V are real and anti-symmetric. 
Define operators 

Q = /- |1)(1| - |2iV)(2iV| and G{u) = {iXi{u) - Ho)-^Q. 

Note that G{u) is well-defined because of Eqs. (3.4,3.8). Representing -01 = Qipi + aijl) + a_R|2A^) 
for some aL,aR G C and rewriting Eq. (3.6) as Q{iXi{u) — -ffo)V'i = uQVipi we obtain 

Q^Jj^ = nG(n)yVi (3.9) 
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and 

i^i = {I- uG{u)Vy\aL\l) + aR\2N)). (3.10) 
Taking into account that \\V\\ < e and using Eqs. (3.4,3.8) one gets 

\\uG{u)V\\ < \\G{u)V\\ < \ <T^^V- (3.11) 

\Xi{u)-l\ 1 + e 

Note that r/ < 1 as long as e < 1. With (3.9), this imphes that ||q;l|1) + aR\2N)\\ = - Q)^Ai|| > 
(1 — J?)!]!/"!!! > 0, hence the expression on the rhs. of (3.10) is non-zero and well-defined because 
of (3.11). We can now use the fact that G{u) and V have support only near the main diagonal to 
show that ipi is a sum of two exponentially decaying states localized near the two boundaries of the 
chain. We will need the following simple fact. 

Proposition 1. Let X be any matrix such that \\X\\ < r] for some < r] < 1 and Xp^q = unless 
\p — q\ < r. Then 

|(/-X)-i| < (l-r?)-!??^ foraUp,q. (3.12) 

Proof. Indeed, let no be the largest integer such that rriQ < \p — q\. Then (X")p^g = unless n > hq. 
It follows that 



i(/-x),]i < ^ i(x«)p,,i < ^ iixr < ^ 



r] 

n=no n=no ' 



□ 



Represent ipi as 

V'l = /^lV'^ + + = 1, (3.13) 
where ip^ and are normalized states defined as 

iP^ {I -uG{u)V)-^\l), {I -uG{u)V)-^\2N), (V^^j-i/;^) = (^-f^l^^-f^) = 1. (3.14) 

(All states and coefficients above depend on u.) Proposition 1 and Eq. (3.11) imply that and 
decay exponentially away from the left and right boundary, respectively: 

\{M'')\<{1-Vr'v'^ and Kp|V^)| <(l-r/)-ir/^ . (3.15) 

Here we used the fact that the normalizing factors in (3.14) are upper bounded by 1. This follows 
using (l|G(n) = which implies ||(/ - uGiu)V)-'^ |l)f > |(1|(/ - uG{u)V)-'^\l)\'^ = 1 and similarly 
for {I-uG{u)V)-^ \2N). 

We would like to show that the states tp^, tp^ can be made real by choosing their overall phase 
(up to exponentially small corrections). Since many equations below include equalities that hold with 
exponentially small corrections, let us set up a special notation. Given numbers x,y G C, we will 
say that x exp-equals y and write x *=? y iff |x — y| < rj^/'^^' for some constant 1 < 7 < 2 and for 
all sufficiently large N. The same notation X Y will be applied to operators X, Y acting on 
meanmg that \\X - Y\\ < ?7^/'>^ for some constant 1 < 7 < 2 and for all sufficiently large N. 

Define a projector Pq = Po{u) as 

^'o = |V'i)(V'i| + lV^i>(V^i|- (3.16) 

The spectral bounds Eq. (3.8) imply that Po{u) has rank 2 and depends smoothly on u. In addition, 
Eq. (3.15) implies that matrix elements of Pq decay exponentially as one deviates from the top-left and 
the bottom-right corner. Let 11^ and 11^ be diagonal projectors that project onto subspaces spanned 
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by the first and the last N basis vectors respectively (the left and the right halves of the chain). 
By definition, + 11^ = /. Define also Pq^ = n^PoII^ and = n^PoII^ (both operators depend 
on u). Using Eq. (3.15) one easily gets 

n^Pon^ =P 0. (3.17) 

This implies that 

Po'= ^'(f + ^(f- (3.18) 

Taking into account that — ^^|| ^ — Y\\ for any operators X, Y whose norm in bounded 
by 1, we conclude that 

{P,^f'='P,^ and {P,y='P,^ (3.19) 

that is, Pq and Pq' are exponentially close to projectors. Furthermore, Po'(O) and Po^(O) are rank- 
1 projectors onto site 1 and 2N, respectively. Continuity then implies that Pq{u) and Pq'{u) are 
exponentially close to rank-1 projectors for all u G [0,u*]. On the other hand, using Eqs. (3.13,3.16) 
and the exponential decay of ip^, if)^, see Eq. (3.15) one gets 

Pr=m\\^''){^''\ + \i>''){i^''\) and Po"" ='l/5i?l'(lV''')(V'''l + lV;'')(^''l). (3.20) 
This is possible only if 

mi^^r^ mi.^)\=' I, and m'^^wr^^. (3.21) 

Let us choose the overall phase oiijj^ and ijj^ to make them real (up to exponentially small corrections), 
that is, {i>^\il)^) =P 1 and (V^^IV'^) =p 1. 

Computing i\i{u) as an expectation value {iI)i\Hq + uV\iIji) and taking into account that i/j^^ip^ 
are (approximately) real, while Hq + uV is a real anti-symmetric matrix, one arrives at 

iXi "= h^RHj^-lHo + n^l^^) + PMij'^lHo + uVli;"^) 0, (3.22) 

since the overlap between ip^ and ip^ is exponentially small, see Eq. (3.15). Recall that this bound 
holds for all u G [0,^*]. Suppose now that u* < 1. Then 

Ai(n*) = (3.23) 

see Eq. (3.7). But this would contradict to Eq. (3.22) for sufficiently large N. Hence u* = 1, that is, 
Eq. (3.22) holds for u = 1. 

The three lowest eigenvalues of Hq + V can be expressed in terms of the Williamson eigenvalues 
oi H + V as 

eI - eI = Xi{l) ^= and - Pj = A2(l). (3.24) 
Weyl's inequality implies that A2(l) > 1 — e, see Eq. (3.5). Hence 

eI-EI = A2(1) - Ai(1) > 1 - e - r?^/2r 

This completes the proof of Theorem 1. 

Finally, let us remark that the real normalized vectors ip^ and ip^ determine zero-energy unpaired 
Majorana modes and localized on the left and the right boundaries by 

2N 2N 

c^='EW)^P and c^='Y.(P\^'')^P^ 
that is, c^c^ + c^c^ 'i? 0, [c^, i^o + V^] '= 0, and [c^, Hq + V] '= (cf. (3.2)). 
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4 Quasi-adiabatic continuation 



In this section we prove Theorem 3 and Corollaries 1,2. For any u E [0, 1] let < \i{u) < . . . < A7v('u) 
be the Williamson eigenvalues of Hq + uV, see Fact 1. The corresponding normalized eigenvectors 
satisfy {HQ+uV)ijjj = iXj{u)il>j. Theorem 1 implies that Xi{u) is exponentially small while Aj(it) > 1— e 
for j > 2 with exponentially small corrections. Define projectors 

N N 

P+ = 5^|V',>(Vil, P-=P+ = ^|V^,)(^,| and Po = \iJi){iJi\ + \i^i){H- 

3=2 j=2 

Here all states and operators are ti-dependent. Fact 1 implies that these projectors form an orthogonal 
decomposition of the identity, Pq + -P+ + -P- = I- 

Proposition 2. Let \ga{u)) be the ground state of Hq + uV in the sector with fermionic parity a S 
{0, 1}. Suppose an orthogonal matrix R S SO{2N) satisfies RP+{0)R^ = P+{u) for some u S [0, 1]. 
Consider any unitary operator U such that 

2N 

UCpU^ = y^^Rp,qCq. 
q=l 

Then U \ga{u)) = |5o-(0)) for o" = 0, 1 up to an overall phase. 

Proof. Let {aj = aj{u)}jLi be the complex fermionic annihilation operators associated with the 

eigenmodes oi H = Hq + uV ., see (3.1). It follows that \ga{u)) can be defined as the state annihilated 
by all dj(u) with j = 2,...,N and having fermionic parity a. The condition RP^[0)R^ = P+{u) 
implies that R'^ipj{u) is a linear combination of ^fc(O) with k = 2, . . . , N for any j > 2. Equivalently, 
is a linear combination of 0^(0) with k = 2, . . . , N. Thus W\gcr{0)) is annihilated by any 
operator aj{u) with j = 2, . . . , N . Finally, the condition det R = 1 implies that U preserves fermionic 
parity, that is, W\g(j{0)) = \g(j{u)), up to an overall phase factor. □ 

We will construct a rotation R = R{u) as above by setting up a differential equation on This 
can be done using standard first-order degenerate perturbation theory. Namely, suppose P+{u + 6u) = 
e^'^"P+(n)e-^^" for some anti-hermitian generator D = D[u) and infinitesimally small 5u. Then 
e~^^"(i?o + uV + 5uV)e^^'^ must commute with P+{u) in the first order in 5u. Hence 

^-^ = [D{ulP^{u)] (4.1) 

whenever D{u) is an anti-hermitian operator satisfying 

Pi{[HQ + uV, D{u)] + V)P+ = 0. (4.2) 

In addition, we would like to find a real solution D(u), since we need to generate an orthogonal 
rotation. Let us choose D{u) as 

D{u)=£{0{V)), (4.3) 

where O and £ are u-dependent super-operators discarding the block-diagonal while keeping the 
'block-off-diagonal' part and inserting the energy denominator. More precisely, 

0{X) = X- Yl Pc^XPa. 
«e{o,+ ,-} 
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To define £ let {It/iA)} be the orthonormal basis of eigenvectors of Hq + uV such that {Hq + uV)\(f)\) = 
\\4>\)- Note that ah eigenvalues A are imaginary. We define £{X) by its non-zero matrix elements in 
this basis, 

{<^mX)\M = -^^^ forA^y. (4.4) 

It is clear that D{u) defined by Eq. (4.3) satisfies Eq. (4.2). It remains to check that D{u) is real and 
anti-symmetric in the standard basis. 

Proposition 3. Let X be any operator on C^^. Then 

0{X)^ = 0{X^), '0{X) = 0{X), (4.5) 

and 

£{X)^=£{X^), I{X)=£(X). (4.6) 

Proof. Identities Eq. (4.5) follow from the fact that Pa are Hermitian and from the complex conjugation 
rules P+ = P_, Pq = Pq. To prove Eq. (4.6) we note that = \(t>-\)- Taking the complex 
conjugate of Eq. (4.4), replacing A and A' by —A and —A' respectively, and taking into account that 
A = —A, A' = —A', one gets £{X) = £{X). The identity £{X)'^ = £{X'^) is obvious since {|0a)} is an 
orthonormal basis and because the eigenvalues {A} are imaginary. □ 

Since V is anti-hermitian and real. Proposition 3 implies that D{u) is also anti-hermitian and 
real. Hence D[u) is a real anti-symmetric matrix. Integrating Eq. (4.1) over u G [0,1] we conclude 
that P+(l) = RP-\-{0)R^ , where R G SO{2N) is the time-ordered exponential describing evolution 
under D{u), 

R = T-exp( [ duD{u)] . (4.7) 



Define an adiabatic evolution Hamiltonian as 

2N 



p,q=l 



Simple algebra shows that [iD{u),Cp] = Yl'^=i ^p,gi'^)'^g- Hence we can choose a unitary U satisfy- 
ing the conditions of Proposition 2 as the time-ordered exponential describing the unitary evolution 
under D(u), namely. 



U = T ■ exp i / duD{u) 
I Jo 

Let us now bound the norm of D{u). We note that 

/>oo 

Jo 

with a similar representation for other blocks PaD{u)Pp with a ^ j3. This yields 
\\P+D{u)Pq\\ < \\V\\ / e-^2Wt+Ai{«)t < \ < 



X2{u) - Ai(n) - 1 - 2e ' 

since \2{u) > 1 — e and \i{u) < e. (In fact, since Ai(n) is exponentially small in N, we could get a 
stronger bound with 1 — e in the denominator which, however, would require to be sufficiently large.) 
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We note that P-D{u)Pq is the complex conjugate of Pj^D{u)Pq and thus ||P_Z)(n)Po|| ^ ^/(l ~ 2e). 
Similar arguments yield ||P+D(u)P_|| < e/2(l — e). Since PaD{u)Pa = for all a, we conclude that 

\\D{u)\\ < \\P+D{u)Po\\ + \\P^D{u)Po\\ + \\P+D{u)P^\\ < . 



This proves Eq. (2.26). 

It remains to show that ofF-diagonal matrix elements of D{u) decay sufficiently fast. We will 
represent D(u) using contour integrals. Let C± and Cq be circles of radius 3e/2 centered at z = ±i 
and z = respectively. Note that C± encircles all eigenvalues ziziXj{u) with j = 2, . . . ,N, while Cq 
encircles ±iXi{u). It follows that 

Pc.{u) = ^f dzizi - Ho - uV)-\ «£{+,-, 0}. (4.8) 

Define Gq{z) = {zl — Hq)^^. We note that ||Go(2:)|| < 2/(3e) for all z G Ca since the spectrum of Hq 
consists of eigenvalues 0, ±1. Hence 

{zl - Ho-uVy^ = {I -uGo{z)V)-^Go{z), where \\uGo{z)V\\<'^ for all z E . 

Furthermore, {Go{z)V)p^q = unless |p — (/I > r + 1 and Go{z)p^q = unless |p — g| < 1. Proposition 1 
from Section 3 now implies that 

\{zI-Ho-uV)jl\<3{2/:i)'-^\\Go{z)\\ < -(2/3)^ 

for all z G Ga- Block-off-diagonality of D{u) implies that D{u) = Ylaj^^ji PaDiu)P/3- Using the contour 
integrals Eq. (4.8) and taking into account Eq. (4.4) we arrive at 

^ f f rl7 rl7' 

P^D{u)Pp = / / {zl -Ho- uV)-'V{z'l -Ho- uV)-^ . 

{2mY J^^c^ J^,(,Cfi z' -z 

Taking into account that \z — z'\ > 1 — 3e > 1/4 we get 

|(P.Z?(n)P;3),,fe| <4(|y 0)' I^P,.l(2/3)^"'^ . 
The bound |i — p| + |^ — fc| > \j — k\ — \p — q\ > \ j — k\ — r yields 



2N oo 

\VpJ{2/3) ''~'''t:'-'' < 2e{2l?,)\^-^\l^'-G\ G = 2^^(2/3 

p,q=l p=0 

One can easily check that (2/3)^^^^ < 1 — l/(16r) for all r > 1. This implies G < 32r. Combining the 
above bounds together we arrive at 

\{P^D{u)Pp)j^k\ < 81 • 2^2 . er2(2/3)l^'-^l/^'^ . 

Combining the six terms PaD{u)Pp together yields the desired bound Eq. (2.27), completing the proof 
of Theorem 3. 

We end this section by proving Corollaries 1 and 2. Using Theorem 3 we can express the perturbed 
ground state in each sector as l^o-) = ^(1)^150-) where U{u) satisfies 

^ = ^D{u)U{u) 
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and where D{u) is a Hermitian effective Hamiltonian with strength roughly e, see Eq. (2.26). Let 

f{u) = {g.\U{u)\g^) (4.9) 

such that 1/(1)1 = |(5(t|5(t)| and let Q = / — \ga){ga\- Taking the derivative over u gives 

df(u) 

= i{g^\D{u)\g,)f{u) + i{g„\D{u)QU{u)\g,) . 

The first term in the derivative only contributes to the phase of f{u) and can be ignored. The second 
term can be bounded using the Cauchy-Schwarz inequality which yields 

< ^{gAb[u)Qb{u)\g^){g,\WQU\g^) = ^l- \f {u)\^ ^/h{^) , (4.10) 

where h{u) = {ga\D{u)QD{u)\ga) . We can bound h{u) by observing that the states QcpCq\go-) and 
Qcp'Cq'\g^) are orthogonal unless {p,q) is a permutation of {p',q')- This yields 



1 1 

^(^) ^ E (I^m(^)I' + \Dp,g{u)Dq,piu)\) = -Tr DHu)Diu) . 



p,q=l 



The bound (2.26) and the assumption e < 1/8 imply ||D(ii)|| < 4e and thus Ti D"^ {u)D{u) < 
2N\\D{u)\\'^ < 32Ne^. We conclude that h{u) < iNe^. Substituting this into Eq. (4.10) yields 

|^Vl-|/WP<26ViV. (4.11) 



Integrating (4.11) over u G [0, 1] and taking into account that /(O) = 1 gives -^/l — |/(1)P < 2ey/N , 
that is, |/(l)p > 1 — ANe^ which proves Corollary 1. 

It remains to prove Corollary 2. By definition, the storage fidelity is 

F\g^{t) = {g\^Mt)){g{t)\)\9), where \g{t)) = e'^^^^+^^' \g) . 

Choose the perturbed ground states such that {ga\gcj) > and let |^) = (l^o) ® |0_r) + l^i) |1_r))/\/2. 
Simple algebra shows that 



II \g){g\ - \~g){g\ ||i < 2^1 - minK^.ls.)^ < 4eViV . (4.12) 

Consider the time-evolved perturbed ground state 

m) = e^(^°+^)* \g) ~ m ® \Qr) + e*'* \gi) ® \Ir))V2 . 
Define also an auxiliary state 

\Ht)) = i\go) \0r) + e''' \gi) \1r))V2 . 

Hence \h{t)) is a ground state of Hq whose time evolution amounts to the dephasing e"'* between Igo) 
and 1^1 ). Using Eq. (4.12) we get 

II \gmg{t)\ - \h{t)){h{t)\ 111 < II \gmg{t)\ - \gmm\ 111 + II \m){m\ - \h{t)){h{t)\ 111 

= 2|||5)(5l -|5)(5ll|i <8e^/iV . 

Replacing \g{t)){g{t)\ by \h(t)){h{t)\ in the expression for (t) and taking into account that \h{t)){h{t)\ 
is invariant under $ec) we arrive at 

mg){t)-\{g\hm'\<8e^ ■ 
This proves Eq. (2.29) since \{g\h{t))\'^ = cos^ {5t/2). 
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5 Lower bound on the storage fidelity 

In this section, we prove Theorem 2. We first express the time-evolved state \g{t)) = e*^^°^^^*|(7) as 
the result of applying a certain unitary to a related ground state \h{t)). This unitary is composed 
of the quasi-adiabatic continuation operator U of Theorem 3 and a time-evolved version thereof. By 
expanding the time-ordered exponential defining U, we rewrite this unitary in terms of correctable 
and uncorrectable errors. Using the exponential decay (2.27) and the localization assumption (2.22), 
we then show that the weight of uncorrectable errors is exponentially small in the system size. This 
implies the desired result. 

5.1 Time-evolution in the Anderson localization regime 

As we argued in Section 2.3, we can formally compute the storage fidelity by considering encoded 
states \g) = aol^o) + ctildi) even though these states are unphysical. Let U be the quasi-adiabatic 
continuation operator satisfying the conditions of Theorem 3, see Eqs. (2.24,2.25), that is, l^o-) = U\ga), 
where \ga) is the ground state of Ho -\-V in the sector with fermionic parity a. We can rewrite the 
time-evolved state as 

\g(t)) = e'(^«+V')tf/e-i(i/o+V)te*(Ho+V')t^t|^) 

^ ^iiH,+V)t^^-i(Ho+V)tu^^h{t)), (5.1) 
where 5 = — is the exponentially small energy splitting, see Theorem 1, and 

\hit))^ao\go) + e^'''aM 

is a 'dephased version' of the initial state \g). (Here and below we ignore the overall phase of quantum 
states). The operator U describes a unitary evolution under a time-dependent Hamiltonian 

. N 

D{u) = ^ X] Dp^q{u)CpCq , < U < 1 , 

p,q=l 

with exponentially decaying interactions, 

\Dp,q{u)\ < cer^AlP-"! where A = (2/3)^/^'' , (5.2) 

see Theorem 3. Defining 
we have 

\g{t)) = U{t)U^ \h{t)) . (5.3) 

We will use (5.2) to argue U (respectively W) has small weight on uncorrectable errors. We will 
subsequently use the localization condition (2.22) to argue that the same is true for U{t)U^. 

5.2 Uncorrectable errors and their qubit weight 

It will be convenient to describe the error correction algorithm of Section 2.3 by mapping the fermionic 
modes to qubits via the Jordan- Wigner transformation 

C2j-i = Zi - ■ ■ Zj^iXj and C2j = Zi - ■ ■ Zj^iYj , (5.4) 
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where j = 1,...,N and Xj (Zj) stands for the Pauh o"^ (cr^) operator acting on a qubit j. This 
estabhshes a one-to-one correspondence between Majorana monomials and A^-qubit Pauh operators. 
For any Majorana monomial E = Cp^ • • • Cp^, we define the fermionic weight f{E) as the number of 
Majorana modes in E, that is, f{E) = k. We also define the qubit weight C[{E) as the number of 
qubits acted on by the Pauli version of E. Applying the Jordan- Wigner transformation Eq. (5.4) to 

the stabilizers Sj = {—i)c2jC2j+i and elementary errors Ej = (— = {—i)c2j-iC2j one gets 

Sj=XjXj+i, 3 = l,...,N-l and E^ = Zj, j = l,...,N. 

The stabilizer code specified by Si, . . . , S^-i is the repetition code on N qubits. It corrects Pauli 
errors of Z-type on any subset of at most N/2 — 1 qubits. Although Pauli errors of X-type cannot be 
corrected, fermionic superselection rules imply that only an even number of X-type errors can occur 
since any X-type error changes the parity of the total number of fermions. Hence admissible errors 
of X-type are always products of stabilizers Sj and can be ignored. A general Pauli error E can be 
uniquely represented as E = E^E^ , where E^ and E^ are products of X-type and Z-type errors 
respectively. Let us say that E is an uncorrectable error iff E^ acts on at least N/2 qubits. It follows 
that any uncorrectable error has qubit weight at least N/2. 

5.3 Errors caused by the quasi-adiabatic evolution operator 

Let us expand the unitary evolution operator {/ of Theorem 3 in the basis of Majorana monomials: 

U = Y,^eE, u;e = ^Tr {Em) . (5.5) 

E 

Here the sum runs over all 4^ Majorana monomials E. The goal of this section is to bound the 
amplitudes oje oi uncorrectable errors E. 

Lemma 1. Let E be any Majorana monomial with fermionic weight f and qubit weight q. Then 

r.2 



lo^i^l < //^A^e^'^, where y ^ , 

1 — V A 

and where the constants c and 1/2 < A < 1 are defined by (5.2). 

Proof. Indeed, expanding the time-ordered exponential defining U we get 



(5.6) 



\uje\ < r max T^^l ^ , 



4'^n! 

n=0 p q 



where ^ 

^i^lu = • • • Dp„,qn{Un)\ " ^ I {E^ Cp^Cq^ ■ ■ ■ Cp^Cq^ 

In the above equations the sums over p and q run over all n-tuples of integers in the interval [1, 2X], 
while the maximum is taken over all n-tuples u of real numbers in the interval [0, 1]. We claim that 



ri"],- = unless 

p,q,u 



2n>f and ^(3 + |pa - gj) > 2q . (5.7) 



a=l 
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The first condition is obvious: the fermionic weight of Cp^Cgj ■ ■ ■ Cp^Cq^ is at most 2n, while has 
fermionic weight f. To derive the second condition we compute the trace by applying the Jordan- 
Wigner transformation to all operators and noting that CpCg has a qubit weight at most {3 + \p — q\)/2. 
Combining the second condition in Eq. (5.7) with the bound Eq. (5.2) results in 



This implies 



P,Q,U — \ J — \ I 



Taking into account that 



we arrive at 



P q \P6Z J 



\^e\ < 2^ — , X 



Al/2)r 



n=f/2 

Since we assumed 1/2 < A < 1, we have A^'^/^ < 3_ Using e" < ef/^e"-/^ and extending the sum to all 
non-negative n we easily get Eq. (5.6). □ 

We next consider the coefficients ve in the expansion 

Uit) = YueE, ^e = ^ Tr(i?t[/(t)) . 

E 

Lemma 2. Let E he any Majorana monomial with fermionic weight f and qubit weight q. Assume 
that the localization condition (2.22) holds for some constants C, ^ > 0, where C > 1. Then 



where y is defined in Eq. (5.6). 



'gJV/(i6.24)gf/i6 i/q>iV/4 and f > iV/24 

max{Al/^e-l/«}^Cfef/8 ifq> N/4 and f < N/24: 

e^/^ i/q<A^/4 



The above lemma is the only step in the proof that uses the dynamical localization condition of 
Definition 1. In fact, the lemma can be proved using a slightly weaker version of Eq. (2.22). Namely, 
we only need a bound 

^E[\Tt{E^it)E)\]<C^e-T 

for any Majorana monomials E,E\ with fermionic weight m and qubit weights q(-E') > N/4^ q(-Ei) < 
This bound can be easily derived from Eq. (2.22) as explained below. 

Proof. We have 

VE = Y,^E,^TT{E^{t)E^) (5.8) 

El 

where A{t) = Qi-{^o+V)t ^^~i{Ho+v)t jg ^.j^g time-evolved operator A. If Ei has fermionic weight m, i.e., 

Ei = Cp^--- Cp^ = Cp> for some pi < P2 < ■ ■ ■ < Pm , (5.9) 
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then 



E^{t) = cp,{t)...cp^{t) = ^ ^ (-l)%^(,,(t)---cp^(„)(t) 



since the time-evolved operators Cp{t) anticommute. With Eq. (2.21), we get 

^ m 



(5.10) 



a=l 



where the sum is over all m-tuples q= (gi, . . . , qm) of integers in the interval [1, 2A^] and Rp^q = Rp^q{t). 
The latter expression shows that Rp^q is antisymmetric in the components of q. In particular, we can 
restrict the sum to m-tuples with distinct entries, qa / qb for a ^h. We conclude that Ei{t) is a linear 
combination of monomials with the same fermionic weight m as, Ei. Using this fact, we can restrict 
the sum (5.8) to monomials Ei with the same fermionic weight f as E. Applying Lemma 1 then gives 

< J2 e'/'A^^^^^e'^^E[|Tr(ii;i(t)ii;t)| 

Si:f(Ei)=f 

'27V 



< 



, je^/^e^^ max E 

f / Ei:f{Ei)=f 



A^(^^)^|Tr(i?i(t)i?t)| 



Er.f{Ei)=f 



A^(^^)^|Tr(i?i(t)i?t)| 



The third inequality of the Lemma follows immediately. Similarly, the first inequality follows because 
gf/8 < gf/i6gAf/{i6-24) foj. f > ^/24. It remains to show that for all Ei with f{Ei) = f < iV/24, we have 



E 



Aq(^i)^|Tr(i?i(t)i?t)| 



< max{A^/^ e-i/«}^C7f if q > iV/4 . 



(5.11) 



As this statement is trivial if q(-Ei) > it suffices to consider the case where c\{Ei) < N/8. If 

El = ioT p = (pi < ■ ■ ■ < pf), we have 



= ^ J^deti?[p,g]c,- 



according to (5.10) and hence 



max E 

^^■q(Ei)<N/8 



X^i^^)l^\Tr{Ei{t)E^)\ 



< max E[|deti2[p,g1|l 

p:q(cff)<N/8 



(5.12) 



assuming that E = c^. Consider two Majorana monomials and with fermionic weight f. We 
claim that 



^\Pa-qa\> N/8 if q{c^) > N/4, q{cjj) < N/8 and f < iV/24 . (5.13) 

a=l 

Indeed, since = (hLi <^Pa) nLi(cpaCgJ (up to a sign), we have 

f 

q(cg) < q(cj?) + J]]q(cp,CqJ . 



a=l 
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Using q(cpCg) < (3 + ||? — g|)/2, we conclude that 

f 



Y,\Pa-qa\> 2(q(s-) - q(cp^)) - 3f , (5.14) 



a=l 



and the claim (5.13) follows. Combining (5.12) and (5.13) with our localization assumption (2.22) 
gives (5.11), as desired. □ 



In summary, Lemma 1 and Lemma 2 show that 

\u:e\ < /3lAf , E[\ue\] < PW2 

where 



if q > iV/4 
1 otherwise 



/3i = A = (2/3)V(4r) Ai = e?^ 

/32 = CeVie a = maxjA^/^, e"!/?, eV(i6-24)| ^ (1 + gi/8)2gj/ _ 

Consider two Major ana monomials Ei,E2 with fermionic weights fi and f2 such that q{EiE2) > N/2. 
We then have 

nc^E.i^E^l] < ^1^2 ■ (AiAa)^ • max{A, A}^ (5.15) 
because q(£;i) + q(£;2) > q{EiE2). 



5.4 Lower bound on the storage fidelity 

We can now easily prove Theorem 2. 

Proof. Represent the time-evolved state \g{t)) in Eq. (5.3) as 

\g{t)) = U{t)U^ \h{t)) = Vgood + ^bad, 

where V'good ^-nd V'bad represent the overall contribution of correctable (qubit weight < A^/2) and 
uncorrectable (qubit weight > N/2) errors, respectively: 

V'good= Yl &EE\h{t)) and Vbad= 9EE\h{t)). 

E:q{E)<N/2 E:q{E)>N/2 

Here 

Ei,E2: 
E2Ei=E 

is the coefficient corresponding to the monomial E in the expansion of the operator U{t)W . We obtain 
the following bound on the norm of V'bad: 

l^E2^Ei\ ■ 

E:q{E)>N/2 Ei,E2: 

q{E2E^)>N/2 
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Taking the expectation over disorder realizations, we get with (5.15), 

fl,f2 ^ ^ ^ ^ ^ 

= ((l + /3i)2(l + /32)=^(AiA2))'^max{A,A}^ . 

Defining 

/(e, C, e, r) = e2^(l + eV8)2(i + ei/4)2(i + ^,i/i6)2 max{(2/3)i/(32-) , e-V€, ei/(i6-24)| ^ (g^^g) 

we conclude that E[||?/;;,ad||] < f{e,C,£,,r)^. We will assume that f{e,C,^,r) < 1. Let V'good ~ 
llV'goodll""'^ • V'good, such that ^good is a normalized state. Simple algebra shows that 

IE[|| Kood) - m) II] < 2E[||V^bad||] < 2/(6,C,C,r)^ . 
Furthermore, using the definition of the error correction map, see Eq. (2.15), we get 

^ec(l4ood>«oodl) = |/i(i))(/iWI 

since C{s)QsE \h(t)) ~ \h{t)) for any correctable error E. This allows us to bound the expected 
storage fidelity as 



E[F{t)] > \{h{t)\g)f - 2V2f{e,C,^,rf = ||ao|' + lail^e^^^f - 2^2/(6, C,C, 



where 5 = E[ — is the exponentially small energy splitting of the ground state, see Theorem 1. 
Since 1^) is a normalized state, the first term on the rhs. is lower bounded by cos^{5t/2). 

Collecting all the constants involved in the definition of f{€,C,^,r) and taking into account that 
e < 1/4 one easily gets the bound 



/(e, C, r) < exp \ar^,/e + /SCeTe 



■ max 



■ _J_ -i 1 

e isr , e « , e 28o 



where a, 13 = 0(1) are some constants. Furthermore, Theorem 1 implies 



, 2e \ 2'- N 
6 < < 2^27. 



Therefore cos'^{6t/2) > 1 — t^2 r . This completes the proof of Theorem 2. □ 



6 Localization and the one-particle problem 

In this section we explain how to map the one-particle Hamiltonian Hq + V corresponding to the Ma- 
jorana chain model to the one-dimensional Anderson model with disorder in the hopping amplitudes. 
We give evidence for the exponential localization of eigenvectors in the latter model by computing 
Lyapunov exponents for the one-particle eigenfunctions. Finally we show how to replace the expecta- 
tion value of the determinant in the multi-point dynamical localization condition Eq. (2.22) by more 
standard multi-point correlation functions, see Section 6.3. 
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6.1 Relation to one-particle Anderson models 



Consider a fixed configuration of the cliemical potentials fii, . . . , fiN and define a triple of one-particle 
Hamiltonians acting on C^''^, C^-'^, and respectively: 



Ho + V 



-/XI 
/XI 1 
-1 



-/X2 

^^2 



1 

-1 —fJ-N 
fJ-N 



(6.1) 











-1 

/i2 

/t2 



(6.2) 



-1 








H. 



+ 

























/i4 






/i4 





IJ-N-l IJ'N-i 



fJ-N 



(6.3) 



Note that Hq + y is the one-particle Hamiltonian describing the Majorana chain model with J = 1 
and |A| = w, see Eqs. (2.5,2.6). The Hamiltonian Eq. (6.2) corresponds to an Anderson model with 
off-diagonal disorder. It describes a particle that can hop between adjacent sites, where every second 
hopping amplitude is random and the remaining hopping amplitudes are homogeneous. Hamilto- 
nian (6.3) involves both off-diagonal as well as (correlated) on-site disorder. Because (6.3) avoids the 
bipartite structure inherent in (6.2), this form will be more useful for our purposes. Furthermore, if 
one neglects the diagonal matrix elements in (6.3), one obtains an Anderson model with uncorrelated 
random hopping amplitudes which has been very well studied. In particular, dynamical localization 
of this model has been established^ in [18]. 

Let us now describe relationships between the three Hamiltonians defined above. We begin by 
noting that Hs and i{HQ + V) can be mapped to each other by a diagonal unitary operator: 



2N 



Hs = i U{Ho + V)U\ where U = Y^ i^'^ 

p=i 



(6.4) 



^Strictly speaking, the dynamical localization bound of Ref. [18] applies to any energy interval that does not contain 
the point E = 0. 



30 



Hence these two Hamiltonians are equivalent as far as their locahzation properties are concerned. 

To describe the map from Hq + V to let us identify the 2N Majorana operators ci, . . . ,C2n 
with standard basis vectors of a Hilbert space Ti = (8) such that 



C2p-1 ^ \p) «) |0) 



and 



C2p ^ \p) CS> |1) , 



P 



1, 



.N 



(6.5) 



We then have 



Ho = T<g)\l){0\-T^ ®\0){1\ 



y = 5^(|l)(0|-|0)(l|) 



S 



T=Y^\k){k + l\ 

k=l 

[f^k\k){k\ . 



N 
k=l 



A straightforward calculation gives 

-{Ho + Vf = (Tt + S)iT + S)(E) |0)(0| + (T + S){T^ + S) (g) |1)(1| 
In terms of these operators, is given by 

H+ = {T + S){T^ + S) - I . 



(6.6) 



The spectrum of is therefore formed by eigenvalues Ej = Xj — 1, j = 1, . . . , N . Localization 
properties of eigenvectors of Hq + V can directly be translated to those of if+ and vice versa. Indeed, 
projecting Williamson eigenvectors ipj onto the even sublattice we obtain a complete family of eigen- 
vectors for H^. Conversely, let (j) ^ be a (real) eigenvector of with an eigenvalue E > —1. 
Using the identity S + 1 = || (T^ + S')|0) p it is easy to check that 



\\{T^+s)m\ 



(rt + 5)|0) |o) + 



ID 



(6.7) 



and i/j are orthogonal eigenvectors of Hq + V with eigenvalues ±i\/ E + 1 (cf. Fact 1). Furthermore, 
since T^^ — 5 is local and the encoding (6.5) respects locality (modulo a factor of 2), localization 
properties of \<j)) translate directly into those of \^) as claimed. 



6.2 Lyapunov exponents 

Since there is a local correspondence between eigenvectors of Hq + V and eigenvectors of H^, see 
Eq. (6.7), it suffices to establish localization properties for the latter. Consider an eigenvector 



N 



H+^ = E^l^ , V' = E(-1)"^"I^)- 



n=l 



Using the explicit form of H+, see Eq. (6.3), one gets 



(m? - i?)Vi 

{ls%-E- l)Viv 



^^Ni'N-l ■ 



n 



2,...,N-1, 



This can be rewritten using the transfer matrix representation as 



Mr, 



tpn-1 



where Af„ 



1 



jJ-n 
IJ-n + 1 





(6.8) 
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It follows that 



Mn--- M2M1 







It will be convenient to introduce an auxiliary transfer matrix 



1 



Then 



■Mn- 



1pn+l 



~1 







1 



Ail 



/J-n 



-^J'n 





G SL{2, 



(6.9) 



■01 V/"n+l, 

It is important that the matrix T„ depends only on Since {//n} are i.i.d. random variables, we 
can regard the product T„ • • • T2T1 as a random walk on SL{2, M). For later use let us denote 



T{C) = 



-c 







Let be the real projective plain, i.e., the set of two-dimensional non-zero real vectors where pairs 
of vectors lying on the same line are identified. Any matrix T G 5-L(2,]R) defines a one-to-one map 
T : P-*^ — )• P-*^. We will use the following statement (see e.g., [17]), which is a variant of Fiirstenberg's 
theorem: 



Theorem 4 ([25]). Let p{nj) be the probability density of Hj. Suppose T(C) G S'L(2,R) satisfies 



iog||r(C)||dp(C) <oo . 



(6.10) 



Let Gp be the smallest closed subgroup of SL{2,M) containing all matrices T{() with p{C) > 0. Assume 
that 

(i) Gp is not compact and 

(a) there is no subset L = {^1,^2} CI P^ with one or two elements such that T{vi) G L and T{v2) G L 
for every T £ Gp. 

Then there exists i > that does not depend on disorder such that 



lim -iog||r„---r2Ti|| 



(6.11) 



with probability 1 over disorder realizations. 



To verify that this theorem can be applied, observe that the integrability condition (6.10) is satisfied 
since we assume that p has compact support [p — r], fi + r]] and ||T'(C)|| < 2 max (C, C~^) (C^ ~ -^)C~^)) 
so the integral in Eq. (6.10) converges absolutely. Let x,y be in the support of p and define 



T ^ T{x)-'T{y) 



a 

a — a'^ a^^ 



where a = x/y . 



Because 



a'' - a ^ a 



G Gp for every p > 1 , 



(6.12) 



(6.13) 
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and we can assume without loss of generality that a > 1, we conclude that Gp is not compact. 
Similarly, we obtain from (6.13) 

lim T^(w) = vi, where vi = 

for any f G assuming a > 1 (that is, x > y). Since all powers of T belong to Gp, we conclude that 
if a subset L C satisfying condition (ii) exists then vi G L. On the other hand, interchanging the 
role of X and y in (6.12) gives a < 1 and leads to 



lim TP{v) = V2 where V2 = ( ^ ) 
p— >oo \ 1 / 



for any v Hence if a subset L C P^ satisfying condition (ii) exists then L = {fi, ^2}- Take x ^ 

in the support of p. Then T{x)v2 = —x ( ^ ] ^ L, hence there is no subset L satisfying condition (ii). 



A theorem by Oseledec [46] and Ruelle [49] states that if a sequence of matrices {T„ G SL{2, R)}„>i 
satisfies 

lim -||T„|| = 

n—>-oo fi 

and (6.11) for some constant £ > 0, then there is a basis {f4.,t;_} of such that 

lim - log ||T„ • • • Tiv±\\ = ±e . (6.14) 

The number i will be called the Lyapunov exponent of the eigenvector ip. From Eq. (6.9) and (6.14) 
we infer that i controls the rate of the exponential growth (decay) of the chosen eigenvector, that is, 

£= lim -log|V„M|. (6.15) 

Let i{E) be the Lyapunov exponent considered as a function of the energy E. The main goal of this 
section is to compute i{E) in the limit of weak perturbations, that is, e = fi + ij <^ 1. We will argue 
that 

mm£{E) ~ } (6.16) 
E log (1/e) 

in the limit e — )• 0, if taking the limit one keeps the ratio h/t] between the strengths of the homogeneous 
and disordered part of the perturbation fixed. Assuming that the localization length ^ scales roughly 
as the inverse of the Lyapunov exponent, this suggests that the localization condition required for 
Theorem 2 holds for sufficiently weak perturbations (cf. Section 6.3). 
Introduce variables 

Zn = —-— for n > 2 and zi = 0. 

Vn 

Rewriting the eigenvalue equations for in terms of Zn we obtain 
zi=0, 

Zn+i = ^ — ' n = l,...,iV- 1, (6.17) 

Zn = • 

P'N 
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The Lyapunov exponent is 



e{E) « ^ XI \zj\ for large N . (6.18) 

We will now use Eqs. (6.17,6.18) to argue that the minimal Lyapunov exponent obeys the scaling (6.16). 

Let us first discuss some numerical evidence for (6.16): we compute i{E) numerically for different 
energies E using the recursion relation (6.17). We choose the chemical potentials identically and 

independently distributed according to (2.7), keeping the ratio n/r] fixed while varying the average 
potential fi (this corresponds to the perturbation strength e). Throughout, we set iV = 10^. Fi gure 6.2 
shows that £{E) is minimal around E y?; furthermore, this minimum becomes more pronounced as 
the overall perturbation strength /i is decreased. Figure 2 shows the value of the minimal Lyapunov 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 



E E 

(a) Lyapunov exponent 1{E) as a function of E, for ^ = (b) Lyapunov exponent t.{E) for ^ — 2 

2~'^ and fi/ri G {1.5,2,2.5}. The inset shows the behavior 
around the minimum. 



Figure 1: The Lyapunov exponent i{E) is positive throughout the interval [—1,1], and achieves its 
minimum around E ^ fi'^. 

exponent i^ain = rainE i{E) as a function of the perturbation strength. We indeed observe a linear 
scaling i^i^ ~ iog\/fi as /i 0, confirming (6.16). 

Figure 6.2 suggests studying the behavior of i{E) for small values of E (i.e., E ~ fj?). In the 
following, we analyze the behavior of the exponent i{E = 0) at zero energy as a function of the 
perturbation strength. Unfortunately, these arguments do not seem to extend to the whole range of 
energies of interest. 

Lyapunov exponent at E' = 

Here we give some heuristic arguments showing that with overwhelming probability over disorder 
realizations, the Lyapunov exponent at energy E = scales as 

£(E = 0) as e ^ . 

logl/e 

This is based on the techniques developed by Eggarter et al [22] which relate the Anderson model to 
certain random walks. 
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Fi gure 2: The niininial Lyapunov exponent ^min 

= miniE i{E) as a function of (log 1/fi) , where fi is 
the average chemical potential. In the limit /x — )• 0, this becomes linear, following (6.16). 



For simplicity we will assume that 

e/2< fin<e for all n (6.19) 

where < e ^ 1 controls the perturbation strength. Specializing (6.17) to E = gives the recursion 
relation 

Z„+l = , Z2 = /i2Mi • (6-20) 

A few starting elements in this sequence are 

1 1 — 

1-Ml l-/i^(l-^l ) 

We will prove that for all sufficiently small n the sequence z„ has alternating signs, namely, z„ > 
for even n and z„ < for odd n. Furthermore, we will show that the behavior of x„ = In (z„) for 
even n can be modeled by an unbiased random walk on a line with a potential barrier on the interval 
(a, oo), where a ~ | loge|. The potential is negligible outside of the interval (a, oo) and grows roughly 
as \z — ap for z S (a, cxd). The walk starts from X2 w | loge| near the left boundary of the potential 
barrier. As soon as the random walk reaches values a;„ — | loge|, it jumps back to its starting point 
reaching Xn+i ~ | log e| in a single step. This jump also breaks the alternating sign pattern of z„ such 
that both Zn and Zn+i are positive. The typical number of steps it takes for the walk to diffuse from 
I log e| to — I log e| is n* ~ I log ep. The random walk picture allows one to partition the sequence {zn} 
into cycles of length roughly | log ep. Within each cycle the sequence {zn\ follows the alternating sign 
pattern. A cycle ends whenever two adjacent z^'s have the same sign. We will prove that each cycle 
increases log \ il)n\ by roughly | log e| . Hence the average increase of log \ il)n\ per step is roughly | log el""*^ 
providing the estimate of the Lyapunov exponent. 

Iterating Eq. (6.20) twice one gets 

Zn+2 = ( ^""o^^"" ] ■ Zn- f3n{Zn), (6.22) 



where 



/?n(.n) = . (6.23) 
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It is important that 

/3„(z) w 1 whenever e < |z| < e"^ . (6.24) 

Furthermore, 

< /3„(z) < 1 for all z > /i„ . (6.25) 

Below we consider only even values of n. Introducing variables x„ = In we obtain 

Xn+2 = Xn + 6^+2 + Gn - 2e„+i + (fn{Xn), 

where 

en = ln|^„/e| and ipnix) = log/3„(e^). 
By definition we have |e„| = 0(1). Simple algebra shows that 

n 

Xn = X2-en- 62 + 2^(-l)-5'' Cj + ifn-2{Xn-2) + ^Pn-i{Xn-A) + • • • + ^2{x2) 
i=2 

for any even n. In order to interpret these equations in terms of a random walk, define auxiliary 
variables Un by 

Vn+i = yn + in + ^2n{yn), where yi = X2 and = 2(e2n+2 - e2n+l)- (6.26) 
One can easily check that 

\x2n — yn\ < 0{1) for all n . 

In addition, ^„ are i.i.d. random variables with E(^„) = and variance 0(1) that can be easily 
computed given the distribution of ^„ . Hence yi,y2, ■ ■ ■ ,yn is an unbiased random walk on a line 
where each step has length 0(1). The term ip2n{yn) can be regarded as an external force, or, position- 
dependent bias. This force turns on only then y„ becomes close to lb I log e|, see Eq. (6.24). In both cases 
the force pushes the walker towards negative values of yn, see Eq. (6.25). In addition, the magnitude 
of the force grows approximately linearly with y„ in the interval y„ G (a, oo) with a ~ | loge|. Indeed, 
in this interval one has 

^2n{y) ~ -log(l + ^^ne^) ^ -\y - \ log/i„||. (6.27) 

Such a force is roughly equivalent to an energy barrier that reflects the walk each time it enters the 
interval (o, oo). From Eq. (6.21) we infer that the starting point of the walk obeys 

yi«|loge|. (6.28) 

Consider now what happens when y„ reaches values of order — | loge|. Equivalently, Zn reaches values 
of order e. Assuming that z„ is positive and satisfies z„ < /i„ we get 

zn+i > ^ ~ . (6.29) 

Equivalently, y = logz^ jumps from — | loge| to +| loge| in a single step. Note also that both z„ and 
Zn-\-i are positive which breaks the alternating sign pattern. We can now start the next walk described 
by Eq. (6.26), now using Zn with odd n. The number of steps needed for each of such walks to diffuse 
from +1 log e| to — I log e| is 

* I 1 |2 

n ~ I log e| . 

Since the walk spends most of its time in the region with negligible external force, we can compute 
the increase of log \ ^pn\ per cycle by setting I3n{zn) = 1. Then we have 

■p^ = exp Xj = exp [Xn] ■ 
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Simple algebra shows that for any even n one has 



Xr, 



I 1 1 

Xj = -Xn + -en + ei = -Xn + 0(1) 



j=2 



Hence we have 



X„«(l/2)y„/2« (-1/2)1 log e|. 



We conclude that the logarithm of increases by roughly 

(1/2)1 log e|. 



per cycle, where cycle consists of n* 
roughly 



£{E = 0) 



1 



7 log I ; 



logep steps. This means that the Lyapunov exponent is 

IV'n-l Xn' loge 1 



n 



loge| 



loge| 



6.3 Localization of multipoint correlation functions 



Here we show how to derive the multi-point dynamical localization condition (2.22) that involves the 
determinant det R[p, q\ from the analogous localization condition that involves more standard multi- 
point correlation functions. Recall that \p — q\i = Yl^=i \Pa ~ Qa\-, where p and q are m-tuples of 
integers. 

Lemma 3. Let p,q be m-tuples of distinct integers in the interval [1,2A^]. Suppose there exists 
constants ci,^i such that 

\p-q\i 



E 



.a=l 



Pa,qa I 



< c'{' exp 



(6.30) 



for all N and for all m-tuples p, q such that \p — q\i > A^/8. Then the multi-point dynamical localization 
condition (2.22) holds for some constants C = O(ci^f) and ^ = 0{^i). 

Remark: The arguments used in the proof of the lemma can also be applied to any fixed disorder 
realization. In particular, if all matrix elements of -R = R{t) obey a bound \Rp,q\ < cie"!^"'^!/'^^ 
for some disorder realization, then Eq. (6.30) holds without the expectation value. The lemma then 
implies that the multi-point dynamical localization holds for the chosen disorder realization with some 
constants C = 0{ci^\), ^ = 0(^i). Repeating all the steps used in the proof of Theorem 2 for a fixed 
disorder realization one would get the same bound Eq. (2.23) on the storage fidelity. 

Proof. For any /? > define a partition function of pairings at inverse temperature (3 as 

m 

Zp{p,q) = e-^^'^(^''') where K(p,g) = ^ |Pa - g.(a)l, 
where Sm is the group of permutations of m objects. The bound Eq. (6.30) yields 



E[|deti2[p,g]|] < J]] E 



.a=l 



'PaSaia) 



<cTYe 



\p-o-{<l)\l/(,l 



cTZy^^ip,q) 



(6.31) 



The statement of the lemma now follows from the following fact. 
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Proposition 4. Let < P < 1 be a constant independent of m and N . Let p, q be m-tuples of 
integers in the interval [l,N] such that pi < p2 < • • • < Pm and qi < q2 < ■ ■ ■ < Qm- Suppose that 
\p-q\i= n{N). Then 

Zi3{p,q) < 0{l/pf"^e-^^^'^'> . (6.32) 

Proof. We claim that it suffices to prove Eq. (6.32) for the special case when all integers in p are even 
and all integers in q are odd. Indeed, Zj^{p,q) = ^/3/2(p'i ^O) where p'^ = 2pa and q'^ = 2qa (the latter 
partition function is defined for a system size N' = 2N). Hence, without loss of generality, all integers 
in p and q are even. Furthermore, shifting each integer in q by ±1 we change the energy Eu(p,q) at 
most by m. This changes the partition function at most by a factor e^"^ which can be absorbed into 
the factor 0(l//3)"'. 

In the following, we prove Eq. (6.32) assuming that all pa are even and all qa are odd. Let us fix 
some permutation a S Sm- We will say that an interval [n,n'] is a pairing iS n = pa, n' = qa{a) or 
n' = Pa, n = q^(^a) ^or some 1 < a < m. We can consider o" as a collection of m pairings. For any 
i G [0, A^] let Xi be the number of pairings that contain the interval [i,i + 1]. Obviously, < Xi < m 
and xq = xn =0. 

Proposition 5. Let E = Efj{p, q). The sequence xq, xi, . . . , xat has the following properties: 

Xi-Xi.i G {-1,0, +1}, (6.33) 

N 

Y^x, = E, (6.34) 

i=0 

N 

^^|xj — = 2m. (6.35) 

1=1 



Proof. Let us prove Eq. (6.33). Suppose i is even, so that i ^ q. Hi ^ p then no pairing can terminate 
at i. Hence a pairing contains the interval , i + 1] iff it contains the interval [i — 1, i] , that is, Xi = Xj_i . 
Suppose now that i £ p, that is there is one pairing that terminates at i. If this pairing lies on the 
left of i, then Xi = Xi-i — 1. If this pairing lies on the right of i, then Xi = Xi-i + 1. The case of odd i 
is analogous. Property (6.34) is obvious since E„{p,q) is nothing but the total length of all pairings 
in a. Property (6.35) follows from the fact that Xi ^ Xj_i only if some pairing terminates at i and the 
total number of points where some pairing terminates is 2m. □ 

We can think of x = {x,} as a step function with 2m + 1 steps. Let hj and wj be the height and 
the width of the j-th step. From Proposition 5 we infer that 

2m 2m 

< hj < m, ho = /i2m = 0, hi = it 1, = E, '^^'Wi = (6.36) 

i=0 i=0 

Note that the locations of the steps are fully determined by p and q since Xi_i ^ Xi iS i £ p or i £ q. 
Furthermore, the number of permutations a G Sm that give rise to the same sequence x = {xj} is at 
most 

2m 

'^(^) = max (/ij, /ij_i) . (6.37) 

i=l 

Indeed, each time we switch from one step of x to another we either terminate or create some pairing 
and we have a freedom to choose which pairing to terminate or create without changing the sequence x. 
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The number of possible choices when we switch from the i-th. step to the (i + l)-th step is at most 
max (/ij, Hence we can upper bound the partition function as 

2m 

Zp{p,q) < Yl e"^" Yl Hmcix {hi, h,.i) . 

n>Eo 0</lo,.--j'l2m<"^ *=1 

ho = h2m=0, 

hi=/ii_i±l 
Y.- Wihi=n 

where Eq = \p — q\i = Q{N) is the ground state energy (one can easily check that for ordered p and q 
the minimum of E(j{p,q) is achieved for the trivial permutation a). Taking into account that Wi > 1 
for all i and using the inequality max (/ij, < hi + 1 we get 

2m / ^ 2m \ ^"i / ^ 2m \ 

max {hi, hi^i) < ( ^ max {hi, hi^i) \ < ( ^ "^i + ^ih J 
i=i \ ^ i=i J V 1=1 / 

Using the constraints ^^1=1 ^« ~ Yl'i=i '^i^i — "-i and the inequality (2m)^"* > {2m)\ we arrive at 

^ „ , , (iV + n)2™ 
j|max(/i„/i,_i) <————. 

i=l 

Let A > be a parameter to be chosen later. Then 

J]max(/.„/l,_l) < A-2-^Mil±:^ = ;^-2mgA(7V+n) ^ 
i=l fe=0 

Since the number of sequences ho, . . . , h2m satisfying conditions Eq. (6.36) is at most 2^"^, we obtain 
Zs{p,q) < (2/A)2™e^^ V e-(^-^)" = ^-(P-x)e,+xn 

n>Eo 

Since we are promised that Eq > aN for some constant a = 0(1), we can choose A = a/3/(2 + a) 
which yields 

Zh a) < (1^^ -a/3(2+a)-l7V 

^{P,q) ^ ^ _ ^_2/3(2+a)-l ^ 

This proves Eq. (6.32). □ 

□ 



7 Numerical simulations 

In Section 7.1, we explain how to efficiently sample from the distribution of syndromes. In Section 7.2, 
we show how to compute fidelities efficiently. In combination, this gives a Monte Carlo method for 
estimating storage fidelities. In Section 7.3, we discuss numerical results obtained using this method. 
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7.1 Simulation of the syndrome measurement 

Let t be some fixed time and let -ir{s) be the probability of measuring a syndrome s = (si, . . . , sn-i) 
on the time-evolved state \g{t)). Our first goal is to describe an efficient algorithm that allows one to 
sample s from the distribution vr(s). An important fact is that both time evolution and the stabilizer 
measurements belong to a class of operations known as fermionic linear optics for which efficient 
simulation algorithms have been described by Knill [42] as well as Terhal and DiVincenzo [53]. For 
example, applying the algorithm of [53] to our settings, we can reduce the problem of sampling s 
from 7r(s) to a series of simpler tasks: sample a bit Sj from the conditional distribution of Sj given 
si, . . . , Sj-i, where j = 1, . . . , N — 1. Using the techniques of [53], the conditional probability of, 
say, Sj = can be computed as a ratio of two determinants representing probabilities of outcomes 
si, . . . , and si, . . . , Sj-i. Once this conditional probability is known, the bit Sj can be set 
by tossing a coin with an appropriate bias. Setting the bits of s one by one starting from si, the 
computational cost of generating one full syndrome sample s in this fashion is 0{N^) since this 
involves the computation of 0{N) determinants of matrices of size 0{N), see [53] for details. 

Here we propose a simplified version of this algorithm in which the computational cost of generating 
one full syndrome sample is only 0(A^^). Our algorithm might also be more stable computationally 
since it avoids computing the ratio of probabilities for exponentially unlikely events. 

Lemma 4. Let p be a fermionic Gaussian state of 2N Majorana modes and 7r(s) = Tr (pQ^) be the 
probability of measuring a syndrome s G {0, l}^^-*^ on p, see Section 2.3. There exists an efficient 
algorithm that takes as input the covariance matrix of p and returns a sample s drawn from the 
distribution 7r(s). The algorithm requires roughly arithmetic operations on real numbers and the 
generation of N — 1 random bits. 

Proof. Define a random sequence of fermionic Gaussian states po,pi,. . . ,pn-i such that po = p and 
Pj is obtained from Pj-i by a non-destructive measurement of a syndrome bit Sj. More formally, let 
Ilj{sj) = (l/2)(/ — i{—l)^i C2jC2j+i) be the projector onto the eigen-subspace with a fixed syndrome 
bit Sj and Pj{sj) = Tv{pj-ill.j{sj)) be the corresponding probability. Then we define 

where Sj is a random bit with a probability distribution (pj(0),pj(l)). Let M^-') be the covariance 

matrix of pj, that is, Mp)^ = (— z/2) Tr(pj[cp, Cg]). Note that the matrix M^"-* is the input of the 
algorithm. By definition, we have 

2p,[sj) = l + {-iy^Mi];^f^, . (7.2) 
Choose any 2j + l<p<q< 2N. Using Eq. (7.1) we arrive at 

= —^T^H-i)CpCgiI - ii-ir^C2jC2j+l)pj-l . 
Pj{Sj) 

Since pj-i is a Gaussian state, we can employ Wick's theorem to compute the trace Tv{c2jC2j+iCpCqPj-\). 
After simple algebra one gets 

(i) 

Here it suffices to only compute matrix elements Mp^q with 

p = 1, 2j + 2 < q <2N and 2j + 2 <p < q <2N. 
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Combining Eqs. (7.2,7.3) we obtain inductive rules for computing conditional probabilities pj{sj). 
The overall computation requires roughly N^/6 steps, where each step amounts to computing the 
righthand side of Eq. (7.3) for some fixed p,q,j. This requires roug hly addit ions, multiplications, 
and divisions on complex numbers. □ 

Let us discuss how Lemma 4 can be used in our setting. Let \g{t)) be the final state of the memory 
before error correction. Taking into account superselection rules we get 

Tr{s) = {g{t)\QMt)) = Yl K\^{9At)\Qs\9At)). (7.4) 

o-=0,l 

Recall that \gu) is the ground state of Hq with fermionic parity a, are amplitudes of the initial 
encoded state, and \ga{t)) = e^i^o+V)t \^g^Y \i follows that 7r(s) is a probabilistic mixture of two 
distributions 7ro(s) and 7ri(s), where 7ro-(s) = {ga{t)\Qs\gcT{t)) ■ The time-evolved states \grT{t)) are 
Gaussian and we can efficiently compute their covariance matrices, namely, M'^(t) = R{t)M'^ {0)R{t)^ , 
where R{t) = e'^"o+V)t ^ S0{2N). Here i^o and V are one-particle Hamiltonians corresponding to 
Hq and V, see (3.3), and M'^(O) = A'l"' are anti-symmetric matrices with non-zero entries (above the 
diagonal) 

^iV = (-ir and M^j^2j+i = ^ for j = 1,...N -1 (7.5) 

The matrix exponential R{t) can be computed in time 0{N^) by finding Williamson eigenvectors of 
Hq + V, see Section 3.1. Applying Lemma 4, we can sample s from Tr^is) and thus we can sample s 
from vr(s). 

It is worth mentioning that an efficient simulation of the syndrome measurement alone can be 
used to assess the storage time of the memory by estimating the probability of 'bad syndromes' as 
defined in [14] (in our context, a syndrome s is bad iff it can only be caused by an error of weight 
roughly N/2). In particular, if the overall probability of bad syndromes is exponentially small (in A^) 
for all t G M, then the storage time of the memory grows at least exponentially with N, see [14] for 
details. 



7.2 Monte Carlo computation of the storage fidelity 

Represent the storage fidelity as an expectation value 

i^,)(t)=E(/,) , 



(7.6) 



where s is drawn from the distribution tt{s) = {g{t)\Qs\git)) and fs G [0, 1] is the fidelity between the 
initial encoded state and a normalized error corrected final state for given a syndrome s. Here the 
expectation value is taken only over the distribution of s (fixed disorder configuration). More formally, 

2 



fs 



\{g\C{s)Qs\gm' _ 1 
7r(s) 7r{s) 



(T=0,1 



(7.7) 



Here we omitted Qs since QsC{s)\ga) = C{s)\ga)- The expectation value E(/s) can be estimated 
with precision 5 by the standard Monte Carlo method. It requires 0(1/5^) independent samples of 
s and computation of fs for each of the samples. The remaining step is to show that for any given 
syndrome s we can compute fs in time 0{N^). 

Recall that the Pfaffian of a complex anti-symmetric matrix A of size 2m is defined as 

(^) = T^;^ Yl '^Sn(r) ^r(l),r(2) " " " ^r(2m-l),r(2m) • 



2'"m! 
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The Pfaffian can be easily computed up to an overall sign using the identity Pf (A)'^ = det (A). It is 
also well known that Pf (A) can be computed directly using the analogue of Gaussian elimination in 
time O(m^). We will exploit the following version of Wick's theorem. 

Proposition 6 (Wick's theorem). Consider an ordered list of operators Li, . . . , L2m such that every 
operator Lj is a linear combination of the Majorana operators ci , . . . , C2n with complex coefficients. 
Let ip he any fermionic Gaussian state. Define an antisymmetric 2m x 2m complex matrix A such 
that Aj^k = (V'l-^i-^fclV') forj < k. Then 

{ij\L,L2---L2mm =Pi{A) . (7.8) 

Proof. Indeed, since {tp) is a Gaussian state, we can regard it is as the fermionic vacuum for properly 
defined complex fermionic modes Op, Op, p = 1, . . . ,N. Rewriting Lj as a linear combination of Op, 
Op, we can now employ Wick's theorem in its standard form. □ 

Let ci, . . . , C2N be the canonical modes of the perturbed Hamiltonian H = Hq + V, see (3.2). As 
was mentioned above, finding the canonical modes amounts to finding Williamson eigenvectors of the 
one-particle Hamiltonian Hq + V. This can be done in time 0{N^). Then 

N N 

^iHt ^ j-j- (cos(Ajt/2) - C2j^iC2j sin(Ajt/2)) = C2j-i (cos(Ajt/2)c2j_i - sin(Ajt/2)c2j) 

= L\ - ■ ■ where -^^2j-i = L^j = cos{Xjt/2)c2j~i — sm{Xjt/2)c2j ■ (7.9) 

Suppose the minimum-weight error consistent with the syndrome s is a product of elementary errors 
at locations ji, . . . ,jq. Then 

C(s) = Ei---Ej^= Y[ {-i)c2j^-iC2j,„ 

m=l 

= L'{ . . . L% where L'i^_, = (-i)c2,,-i , L'^^ = C2j, . (7.10) 
With (7.9) and (7.10) we conclude that there are operators Li, . . . , L2m, m = N + q such that 

{9a\C{s)\ga{t)) = {g„\Li ■ ■ ■ L2m\9a) , 

where each operator is a known linear combination of the Majorana operators cj. Expectation values 
-^j,k = {ga\LjLk\ga) can be easily computed using the covariance matrices M"", see Eq. (7.5). Hence 
we can compute {g(^\C{s)\ga{t)) using Proposition 6 in time O(iV^). 

The missing ingredient to compute fs is the probability vr(s), see Eq. (7.7). Taking into account 
Eq. (7.4), it suffices to compute the overlaps {ga{t)\Qs\9a(t)) . Since \ga{t)) has fermionic parity a, we 
conclude that Qslgait)) is an eigenvector of {—i)ciC2N with an eigenvalue (— 1)'^"'"°'(*\ where a{s) is 
number of non-zero syndrome bits modulo two. Hence (5'o-(i)|0s|5o-(i)) = l(5o-(t)|V's)p5 where \ips) is 
a Gaussian state with a covariance matrix such that the only non-zero elements of above the 
diagonal are 

(Mf )2j,2i+i = (-1)'^^- and (Mf )i,27V = (-l)'^+'^(^) . 
Using the standard formula for the overlap between two Gaussian states one arrives at 

{9At)\Qs\9a{t)) = 2-^v'det (M-(t) + Mf) . 

To summarize, we can compute fs using Eq. (7.7) in time 0{N^). 
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7.3 Numerical results 



We have numerically simulated chains of various sizes in the range 4 < < 256 with and without 
disorder. Throughout, we consider the Hamiltonian (2.1) with w = |A| and J = 1. We compute 
the storage fidelity F\g){t) for the maximally entangled state \g) = {\go) ® |0/j) + l^i) (g) |1_r)) and 
corresponding storage times Tstorage(-^o) for different thresholds Fq close to 1. 

For system sizes up to iV < 12, we compute F^g)it) exactly by summing over syndromes, while 
for larger system sizes, we sample syndromes 10^ times for each time t to estimate the fidelity. Due 
to the A^^-scaling of this algorithm, this takes on the order of a day of a single desktop computer 
for iV = 256. 



Clean case 

We first consider the clean case with perturbation given by a uniform chemical potential fi G {0.1, 0.13, 0.16} 
and = 6. As the perturbation strength is ||y|| = e = /i, we have A^ <C l/e^ and we expect the 
storage fidelity to be dictated by dephasing due to the exponentially small energy splitting 6, see 
Corollary 2 in Section 2.4. Indeed, we recover the behavior F^g){t) ~ cos^{6t/2) as shown in Fig. 3(a). 
Fig. 3(b) illustrates how this simple oscillatory behavior disappears for larger e. 




(a) In the regime e^A^ ^ 1, the behavior of the storage 
fidehty F^g^{t) is dictated by dephasing due to the expo- 
nentiaUy smaU splitting of the ground state degeneracy. 
This figure shows the error- corrected fidehty superimposed 
with the function cos^{St/2), where S is the gap of H. The 
system size is A'' = 6. 



^ 0.85 




10 20 30 40 

t 

(b) When e^iV > 1, the storage fidehty f|g)(t) has a more 
complex behavior. Here the system size is N = 12. In this 
and the following figures, we use lines as a guide to the eye. 



Figure 3: (No disorder) The storage fidelity -F'|g)(i) for different perturbation strengths e in the regimes 
e'^N < 1 and e'^N > 1. 

Fig. 4(a) shows how the storage fidelity -F|g^ (t) depends on the system size A^ for a fixed perturbation 
strength (e = 0.7). Figure 4(b) is derived from the data of 4(a) and shows the storage time Tstorage(-Pb) 
in the clean case as a function of log2 A^ (for different fidelity thresholds Fq). We observe a logarithmic 
scaling (cf. (1.4)) in the clean case with a sufficiently strong perturbation. 
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(a) The storage fidelity F^g^ (t) as a function of time, for 
different system sizes but fixed perturbation e = 0.7. 
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(b) Storage times rstorage(fb) for different fidelity thresh- 
olds Fo as a function of log2 A'^. Also shown are linear fits to 
the values at = 16, 32, 64, 128. For smaller system sizes, 
the storage time deviates from its asymptotic logarithmic 
scaling, presumably due to boundary effects. 



Figure 4: (No disorder) Fidelity and storage time vary with system size for a fixed perturbation 
strength e = 0.7. 




Figure 5: (No disorder) Storage times vary with system size, for different perturbation strengths e S 
{0.5,0.7} and fidehty threshold Fq = 0.95. 

This confirms predictions made recently by Kay [36] based on a mean-field analysis of error cor- 
rection. Following [36] we can describe the time evolution of the encoded state \g) = {\go) + \gi))/\/^ 
as a quantum quench in the transverse field Ising model. Indeed, applying the Jordan- Wigner trans- 
formation, see Section 5.2, the perturbed Hamiltonian becomes 

N-l ^ N 
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while the initial encoded state becomes the tensor product of |+) states, 



\9) = 1+)^^. 

Each elementary error Ej = Zj flips the corresponding qubit mapping |+) to |— ). The average fraction 
of flipped qubits in the final state is (1 — m^{t))/2, where m^{t) is the magnetization, 

"^'(*) = N Y.(9it)\XMt))^ m) = e^(^°+^)* \g). (7.11) 
i=i 

The time evolution of m^{t) in the homogeneous case {jjij = fi for all j) has recently been computed by 
Calabrese et al [12] for quantum quenches within the ferromagnetic phase, that is, < 1. In this regime 
m^{t) was shown to decay exponentially with time, m^{t) ~ e~^*, where the coefficient c > can be 
easily computed from the quasiparticle spectrum and Bogoliubov angles, see [12] for details. Within 
the mean-field approximation correlations between errors on different qubits are neglected. Then the 
number of errors in the final state has expected value A''(l — m^{t))/2 and standard deviation 0{-s/N). 
Hence the error correction is likely to succeed whenever 



A^(l - rrf{t))/2 + O(VA^) < N/2, 

that is, 



This gives the desired logarithmic scaling of the storage time. 



Disordered case 

Next we choose site-dependent chemical potentials {/Xj} according to (2.7) with {xj} drawn indepen- 
dently and uniformly on [—1,1]. To compare with the clean case, we use A'^ = 12 and the average 
chemical potential /i = 0.5 to be in the regime e^N > 1 of intermediate system size (in the clean 
case where e = /i), see Section 2.4. Fig. 6(a) shows the behavior of storage fidelity F\g){t) without 
disorder [r] = 0) with disorder strength r] = 0.25 for several disorder realizations. 

Fig. 6(b) shows that the (disorder-averaged) storage time scales roughly linearly with the system 
size in the presence of disorder. This is in stark contrast to the clean case (Fig. 4(b)). While this 
figure clearly shows a significant enhancement of storage times, it falls short of the exponential scaling 
of Theorem 2. It appears that the regime where Theorem 2 applies is outside the range of what is 
numerically accessible. 



Pseudorandom disorder 

Suppose now that disorder represents a controlled external potential introduced on purpose to enhance 
the storage time of the system. Our numerical results demonstrate that for fixed parameters of 
the model the storage time varies strongly for different disorder realizations lacking a self-averaging 
behavior. It shows that some disorder potentials are better than others in terms of their ability to 
suppress propagation of excitations and enhance the storage time. A natural question is whether it 
is possible to choose a deterministic disorder potential that typically outperforms the random one in 
terms of the storage time when combined with an unknown homogeneous perturbation. 

It is well known that certain deterministic potentials can give rise to localization of eigenvectors in 
the Anderson model, mimicking random disorder. Well-studied examples include periodic potentials 
whose period is incommensurate with the one of the lattice, such as Harper's potential [47, 23] and 
pseudo-random potentials [26, 11]. 
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(a) Storage fidelity F\g) (t) for the clean case with jj. — e — 
0.5 and for 5 disorder realizations with /i — 0.5 and ran- 
domness strength -q = 0.25. The system size is = 12. 



X 77 = 0.25 

o average over realizations 




2 4 6 8 

log,iV 

(b) Plot of logj Tstorage(^o) versus log2 A'' for the fidelity 
threshold Fq = 0.96 with the same parameters as in 
Fig. 6(a). For each system size, 10 different disorder re- 
alizations are considered. A straight line is fitted to the 
average logj T over disorder realizations, showing a linear 
relationship between Tstoragc(-fo) and A'^. 



Figure 6: (Disorder) Fidelity and storage time for different randomness realizations and system sizes, 
for fixed perturbation- and randomness strength 



Here we focus on a potential generated by iterations of the logistic map, namely, fj,j = iJ,+r]{l — 2yj), 
where the sequence {yj} is defined by the initial condition yi G [0, 1] and the rule 

yj+i = a-yj{l-yj). (7.12) 

One can easily check that yj G [0, 1] for all j provided that < a < 4. The sequence {yj} is known to 
exhibit chaotic behavior for 3.57 < a < 4 and almost all initial conditions. 

With the hope of getting a long storage time, it is desirable to select values (yi,a) giving rise to 
highly localized wavefunctions. To define an empirical measure of localization, we use the symmetric 
matrix Hg of Eq. (6.1). This matrix has eigenvectors \(l)f) with eigenvalues ±Aj, where Ai < . . . < 
Xn are the Williamson eigenvalues of Hq -t V. Omitting the eigenvectors ) corresponding to 
the boundary modes, we expect the matrix entries of Tp^g = Y2a=2 se{+ {^011)1 to decay 

exponentially away from the diagonal, Tp^q < c- e"'^"'^'/^. Accordingly, we define ^{p) by least-square 
fitting a straight line to (g, logTp^g) such that 

Tp,q ^ c{p) ■ e'lf-^l/^'^P) for all g > iV if p < iV 
Tp,q ^ c{p) ■ e-|P-^l/«'(P) for all g < TV if p > N . 

We then define the effective localization length as ^eff = ^^^p^{p)- This quantity can be computed in 
time 0{N^) from the sequence of on-site potentials {fj.j}. 
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(a) This shows Tstorage(Jb) compared to log^off for ran- (b) Plot of logj T for the average over disorder reahzations 
dom and pseudorandom disorder. The rightmost box cor- and of fogj Tstoragc for certain pseudorandom sequences. 



responds to a choice of pseudorandom disorder parame- 
ters equal to {yi,a) = (0.2845,3.9914) (cf. (7.12)). As 
expected, choosing these parameters so as to minimize the 
effective localization length ^cff appears to be a good strat- 
egy for maximizing storage times. Here A'^ — 64, /i = 0.5, 
?7 = 0.25 and Fq = 0.97. 



Figure 7: (Pseudo-random disorder) Enhancement of storage times over random disorder 

We have computed ^eflf and storage times for N = 128 and a discrete set of initial values {yi,a). 
We expect initial values resulting in small values of ^eff to result in long storage times. Figure 7(a) 
qualitatively confirms this expectation: it shows that ^eff indeed provides a crude measure for the 
quality of the quantum memory. Figure 7(b) shows the behavior of the storage time for two different 
choices of {yi,a) and different system sizes. It shows that these pseudorandom sequences give better 
storage times than the average over truly random realizations. 
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